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A. V. Sokolov, A. A. Krasnov, N. V. Kuz’mina, and Yu. F. Stus’ 


Abstract 


Methods to measure absolute gravity and deflections of the vertical on a moving base are 
presented. The breadboard of integrated gravimetric system is described. The first results 
of experimental studies confirmed the possibility of high-precision measurement of the 
absolute values of gravity and deflection of the vertical at sea. 


Keywords 


Absolute gravity - Deflections of the vertical - Gravity measurements - Integrated gravi- 


metric system 


1 Introduction 


Knowledge of absolute values of gravity and deflection of 
the vertical (DOV) is essential for solving a number of 
problems in geodesy, high-precision inertial navigation and 
fundamental position, navigation and time support. In water 
areas, the absolute values of gravity are calculated as a sum 
of the gravity a priori value measured at an onshore reference 
station, and of the observed gravity measured with relative 
gravimeters from marine and air vehicles. The errors in 
determining the drift and scale of relative gravimeters reduce 
the accuracy of determining the absolute values of gravity, 
and the necessity to regularly reference the results of offshore 
measurements to the onshore absolute value imposes serious 
operational limitations on the gravity survey procedure. In 
view of the fact that the geological exploration tasks require 
only the knowledge of character of gravity anomalies in 
the survey area, most of modern geophysical works are 
carried out without precise referencing of measurements to 
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the reference stations, which considerably increases the error 
of absolute values of gravity measured in water areas. 

There is a known technique for the relative gravimeters’ 
measurements re-calculation to the absolute level, using 
the global geopotential models (Zheleznyak et al. 2015). 
However, limited spatial resolution of the global models, as 
well as their significant errors in the areas with high gravity 
anomalies can make it impossible to precisely determine 
the absolute values of gravity using the above technique. At 
present, there are no commercial devices for measuring the 
absolute value of gravity from moving vehicles. A number 
of companies are currently studying and developing such 
devices (Bidel et al. 2018; Baumann et al. 2012). 

Integrating the gravity anomalies according to Vening- 
Meinesz formulae has been the main method of high- 
precision determining of DOV in water areas for as long as 
100 years (Vening-Meinesz 1928). However, this method is 
extremely labor-consuming, since it requires accumulation of 
background gravimetric data for an area that is much larger 
than the point specified for determining the DOV. Moreover, 
the error of the absolute value of gravity discussed above is 
a methodological error of the DOV calculating according to 
Vening-Meinesz formulae. Also it should be noted that the 
existing methods of space geodesy do not provide the DOV 
precisely enough, especially in regions with large gravity 
anomalies (Koneshov et al. 2013). 

A design concept of an integrated gravimetric system 
for measuring the absolute gravity on a moving base, and 
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the first results of its bench tests were discussed at the 
Symposium TGSMM-2016 (Peshekhonov et al. 2016). Such 
system is useful to improve accuracy of gravity surveys that 
were carried out without reference measurements. Further 
development of the concept resulted in the integrated system 
construction and incorporation in the equipment for the DOV 
direct measurement by astrogeodetic method at sea which is 
essential in inertial navigation (Chatfield 1997). 

A breadboard of gyrostabilized zenith telescope was 
made, and its operation methodology was developed. This 
paper presents the principles of the integrated system 
construction, including an absolute gravimeter, a zenith 
telescope, a system for their gyroscopic stabilization, and 
receiving equipment of global navigation satellite systems 
(GNSS). The integrated system is intended for measuring the 
absolute values of gravity with 1 mGal RMSE and DOV with 
1 arcsec RMSE at sea. The system operation methodology 
is described, and the results of field testing of the system 
breadboard are discussed. 


2 Principles of Measuring the Absolute 
Values of Gravity and DOV at Sea 


Land gravity instruments for high-precision determination of 
the absolute values of gravity are commercially available and 
are based on measuring the time and length intervals of a 
test body fall in vacuum (Vitushkin 2015). Offshore appli- 
cations of such devices are limited by the effect of inertial 
accelerations, caused by pitch/roll and orbital motion, on the 
measuring system of the gravimeter. For this reason, in case 
of moving vehicles, measurements are taken using relative 
gravimeters with a gyroscopic system for the sensitive ele- 
ment stabilization in the horizon plane. At the same time, the 
influence of the vertical component of inertial acceleration is 
compensated by low-frequency filtering methods involving 
the external navigation data (Stepanov and Koshaev 2010). 
Obviously, absolute measurement of gravity at sea 
requires the sensitive axis of the absolute gravimeter to 
be stabilized in the direction of the local vertical with an 
accuracy of about 15 arcsec, in order to remove the effect of 
horizontal inertial accelerations. However, due to the effect 
of vertical inertial acceleration on the absolute gravimeter, 
the vast majority of its measurements are unreliable, and the 
methods of frequency filtering cannot be used in absence 
of a model of the test body motion in the field of inertial 
accelerations. The solution to this problem was found by 
integration of the initial data of the absolute and relative 
gravimeters. The idea of data integration is as follows. Based 
on the relative gravimeter data, current inertial vertical 
accelerations are determined, and the absolute gravimeter 
measurements at which the accelerations were minimal are 
selected. The resulting measurements of the absolute value 
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of gravity are used for compensating for the errors in the 
relative gravimeter, caused by drift and nonlinearity of the 
scale, as well as the measurements referencing to the absolute 
level (Peshekhonov et al. 2016). 

Field digital zenith cameras have been developed and 
widely used for measuring another important geodetic 
parameter DOV on the ground (Hirt and Biirki 2002; Tian 
et al. 2014; Halicioglu et al. 2012; Gerstbach and Pichler 
2003). They help to implement the astrogeodetic method 
of DOV absolute measurement based on determining the 
direction to the stars located at the zenith, with known 
equatorial coordinates (right ascension a and declination 8). 
At that, the equivalence of astronomical coordinates (latitude 
~, longitude à) of the observation point and the equatorial 
coordinates of the stars is used. This equivalence is due to 
the validity of the following relations (Torge 2001): 


where 9 is Greenwich apparent sidereal time.The Helmert, or 
astronomic, DOV are calculated in accordance with the basic 
expressions (Jekeli 1999): 


where & is DOV projection on the meridian plane; n is DOV 
projection on the prime vertical plane; B, L are geodetic coor- 
dinates of the observation point (latitude and longitude).As 
with the absolute gravimeter measuring absolute values of 
gravity at sea, the zenith telescope obviously needs to be sta- 
bilized for measuring DOV on a moving vehicle. However, 
to ensure the DOV measurement accuracy at the level of few 
seconds of arc, the required precision of stabilization should 
be comparable. Even on a slightly oscillating base, it is 
almost impossible to achieve such precision of stabilization 
by gyroscopic means only. Therefore, the only solution is to 
stabilize the telescope’s sight axis along the normal to the 
reference ellipsoid rather than in the direction of the local 
vertical. The above normal is drawn according to the data 
from two sources: the zenith telescope itself, and the GNSS 
receiver (Peshekhonov et al. 1995). 

A schematic of a gyrostabilized zenith telescope is shown 
in Fig. 1. Stabilization motors are controlled by the sig- 
nals of mismatch between the geodetic coordinates and the 
astronomical coordinates of the crossing point of the zenith 
telescope’s sight axis with the celestial sphere. 

A block diagram of the algorithm for determining the 
DOV components on a moving base is presented in Fig. 2. 
The objective of near-zenith stars observation is to record a 
sequence of frames containing images of the stars, using a 
TV camera, with the frame record time being simultaneously 
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Fig. 1 Schematic of a 
gyrostabilized zenith telescope 


Stabilization /] 
motor 


fixed. In each frame, the coordinates of energy centers of the 
stellar images are determined and then used, along with the 
star catalog data, for stars identification (Mantsvetov et al. 
2006). 

As a result of identification, a data set is formed, where 
the coordinates of stars on the images are matched to the 
equatorial coordinates of stars from the catalogue. Based on 
this set, the parameters of transformation between the frame 
of the TV camera photodetector and the standard frame are 
calculated; these parameters are used for transforming the 
coordinates of the photodetector central point into equatorial 
coordinates and then, taking into account Greenwich appar- 
ent sidereal time, into astronomical coordinates. Moreover, 
the transformation parameters are useful in determining the 
angle of the photodetector frame rotation relative to the stan- 
dard frame (azimuth of the TV camera row), which is used 
for feedback in the azimuth stabilization loop. After that, 
deviation of the astronomical coordinates of the crossing 
point of the zenith telescope’s axis of sight with the celestial 
sphere relative to the geodetic vertical is determined. The 
deviation values are taken into account in the readings of 
accelerometers, and also used in the stabilization loops. Then 
linear accelerations are subtracted from the accelerometer 
readings using GNSS receiver data. 

In order to compensate for the tilt of the sight axis relative 
to the rotation axis of the optronic device, and to prevent 
the effect of accelerometer bias, observations are carried 
out in two diametrically opposed positions, with the zenith 
telescope being turned to 180 degrees. The final values 
of the DOV components are determined by averaging the 
accelerometer data obtained in two positions of the zenith 
telescope. 


Zenith 
telescope 


GNSS 
receiver 


Stabilization 
motor 


3 Integrated System Operation 
Methodology 


Combined operation of an absolute gravimeter and a relative 
one does not suggest any serious changes in the methodology 
of marine gravity surveys. The gravimeters are placed on a 
common base onboard a vessel, as close to the vessel’s meta- 
center as possible (Fig. 3). The navigation data, including 
geographic coordinates, time stamps, and the data from the 
absolute and relative gravimeters are received to a common 
laptop which is also used for further post-processing of the 
data. 

The survey starts from taking reference measurements at 
port, to determine the absolute value of gravity according 
to the methodology described above, as well as the initial 
value of the drift of the relative gravimeter. The gravity incre- 
ments in the survey profiles are measured with the relative 
gravimeter. If the sea is weak (sea state 0,1,2), it is possible 
to determine the absolute value of gravity; to do so, the vessel 
should stay at the specified point of the water area for one to 
two hours. The survey results are processed in offline mode. 
Some additional offshore reference points would provide 
more detailed estimate of the relative gravimeter’s drift and 
reduce the effect of the scale factor error during operations at 
a significant distance from the reference gravity station. 

The zenith telescope is installed on the gyrostabilizer 
instead of an absolute ballistic gravimeter. In view of the 
fact that DOV measurements are taken at night and only 
in case weak sea and clear sky, this circumstance does not 
impose any significant limitations on the methodology of the 
integrated system operation during gravity measurements. 


Accelerometer Stellar sky 
readings image 


Time of image 
registration 


Identification of stars 


Determination of parameters for transforming the coordinates 
of stars inimage to standard coordinates 


Determination of the astronomical coordinates of the point 
corresponding to the intersection of the AZT rotation axis with 
the celestial sphere 


Correction of astronomical coordinates for Earth polar motion 


Calculation of mismatch error of ZT sight axis and geodetic 
zenith 


Correction of accelerometer readings 
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Two diametrically opposed positions of Zenith Telescope 


Star Geodetic 
catalogue coordinates 


Converting the accelerometer readings 


Averaging the accelerometer readings and calculating DOV 


Fig. 2 Block diagram of the algorithm for determining the DOV components on a moving base 


4 Integrated System Breadboard 


To verify the proposed methods for measuring the absolute 
values of the gravitational field parameters, a breadboard of 
the integrated system was assembled and tested. The bread- 
board composition was based on commercially available 
equipment. 

For continuous measurement of gravity increments, 
the breadboard included a mobile gravimeter Chekan-AM 
(Shelf-E model) designed by Concern CSRI Elektropribor, 
JSC (Krasnov et al. 2014; Evstifeev et al. 2014). This 
gravimeter is based on a gravity sensor with double quartz 


elastic system, installed in a two-axis gyrostabilizer. Root- 
mean-square (RMS) error of gravity increment measurement 
with the gravimeter Chekan-AM under the action of dynamic 
disturbances does not exceed 0.4 mGal (Sokolov et al. 2016). 

In the breadboard, the absolute value of gravity is mea- 
sured with a ballistic gravimeter GABL-PM designed for 
field operation by the Institute of Automation and Elec- 
trometry of the Siberian Branch of the RAS (Bunin et al. 
2010). RMS error of measurement of the gravity with the 
gravimeter does not exceed 5 Gal. The dimensions of the 
gravimeter GABL-PM is 42 x 47 x 93 cm and its weight is 
59 kg, which made it possible to install it in the gyrostabilizer 
commercially manufactured by Concern CSRI Elektropribor, 
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Fig. 3 A schematic of integrated system 


JSC. The operating principle and the sensitive elements of 
this gyrostabilizer are similar to the stabilization system of 
the gravimeter Chekan-AM, while its weight and size are 
several times greater, so it could be used for the absolute 
gravimeter stabilization. 

The main elements of the zenith telescope are a catadiop- 
tric lens with a focal distance of 2 m, and a specially designed 
TV camera KT-62 with intrinsic function of automated 
determination of the coordinates of energy centers on the 
images of point objects. Such a system of a lens and a 
camera can determine the position of point objects on the 
TV camera plane with an accuracy of 1/20 pixel, taking 
into account the influence of various external factors, which 
corresponds to 0.1 arcsec. To determine DOV, electronic 
inclination sensors Zerotronic (accelerometers) with an error 
of 0.2 arcsec are installed on the zenith telescope base. The 
error in determination DOV components for such system is 
less than 0.5/. The zenith telescope has an intrinsic azimuth 
rotation drive, and due to its configuration it can be installed 
in the gyrostabilizer instead of the absolute gravimeter. 

To determine the geographic coordinates of location and 
to receive the time stamps, the breadboard included a GNSS 
receiver Javad of geodetic accuracy grade. 


5 Results of the Bench Tests 
Bench tests of the gravimetric system breadboard were car- 
ried out at the Concern CSRI Elektropribor. The test program 
included tests on the breadboard system accuracy on a fixed 
base as well as on dynamic benches of vertical displacements 
and a sea wave simulator. 

The tests of the breadboard instrumental accuracy were 
carried out on a fixed base during 3 days. The methodology 


of continuous measurement of the absolute gravity consisted 
in a series of initial measurements by GABL-PM, continuous 
measurements of the gravity increments by Chekan-AM, and 
a series of final measurements by GABL-PM, which were 
used to refine the drift of the Chekan-AM gravimeter. The 
measurement results are shown in Fig. 4. 

The RMS error of the initial and final measurements taken 
with the GABL-PM gravimeter was © gabs = 0.02 mGal. After 
taking into account the drift of the Chekan-AM gravime- 
ter using the GABL-PM gravimeter readings and after the 
introduction of correction for the lunar-solar tide effects, 
calculated theoretically, the RMS error of the Chekan-AM 
gravimeter measurement was Ogre) = 0.03 mGal. Since the 
measurements of the gravimeters GABL-PM and Chekan- 
AM are independent, the instrumental accuracy of contin- 
uous measurement of the absolute gravity obtained with 
the mock-up integrated system can be estimated from the 


formula: 
Ogm = YOZ + FF.) = 0.04mGal. 


During three nights we put the zenith telescope on the 
gyrostabilizer and provided measurements outside the test 
room. The results of DOV measurements are shown in Fig. 
5. 

Standard deviation of the DOV measurements did not 
exceed 0.57. So installation of the zenith telescope on the 
gyrostabilizer do not affect on the accuracy of the system on 
the fixed base. 

Another aim of the bench tests was to estimate the 
accuracy of absolute gravity measurements on a weakly 
oscillating base. We didn’t provide dynamic tests of the 
gyrostabilized zenith telescope because it is impossible to 
observe stars at the test room. 

As can be seen from Fig. 6, the spread in the GABL-PM 
readings is 200 mGal even with minimum heaving with an 
amplitude of 0.2 m and a period of 200 s, which additionally 
confirms the necessity of using both types of gravimeters for 
occasional determination of absolute gravity aboard a vessel. 

The data from the relative gravimeter were used to deter- 
mine vertical accelerations Wz, which made it possible to 
choose from a set of discrete gravity measurements of the 
absolute gravimeter, the measurements that were performed 
under the best conditions. The criterion for choosing g 
measurements was the limiting value of vertical accelerations 
of less than 3 mGal. It took about 2 h to obtain a set of 100 
reliable measurements of g at heaving with an amplitude of 
0.2 m and a period of 200 s, and the standard deviation of the 
GABL-PM gravimeter readings was 0.91 mGal. However, 
as can be seen from Fig. 7, the difference in the mean 
values of the measurements of the GABL-PM gravimeter 
before, during and after heaving by modulus was 0.11 mGal, 
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Table 1 Estimation of the breadboard gravimetric system accuracy 


| Bench of vertical displacements 
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| Sea-wave simulator 


Test series | Fixed baseg (o,,), mGal 
No. 1 | 982000.01 (0.04) 
No. 2 | 982000.03 (0.04) 


Staticsg (0), mGal 
982000.02 (0.05) 
982000.01 (0.05) 


| Dynamicsg (0), mGal 
| 982000.13 (0.91) 
| 982000.09 (0.88) 


| Staticsg (0), mGal 
| 982000.03 (0.05) 
| 982000.02 (0.05) 


| Dynamicsg (og), mGal 
| 982000.11 (0.75) 
| 982000. 12 (0.78) 


Fig. 8 System breadboard aboard the vessel 


which was the first reliable demonstration of the feasibility of 
measuring the absolute gravity on a weakly oscillating base. 

The final stage of the bench tests was the breadboard 
testing on the sea-wave simulator, which allows specifying 
three-component angular motions of a mobile platform with 
the equipment under test fixed to it. The method for estimat- 
ing the breadboard accuracy at the sea-wave simulator was 
identical to that described above for testing at the bench of 
vertical displacements. The results of two series of mock-up 
tests on dynamic benches are summarized in Table 1. 

From the results of the bench tests it can be seen that 
imitation of smooth sea waves does not significantly affect 
the measurements of the absolute gravity with the breadboard 
gravimetric system. The standard deviation of the heaving 
measurements was less than 1 mGal, and the difference in 
the mean values in statics and dynamics did not exceed 
0.1 mGal. Based on the positive results of the bench tests, we 
proceeded to the next stage — the sea trials of the breadboard 
gravimetric system. 


6 Results of the Sea Trials 


The system breadboard was tested at Priozersk seaport and 
in the Ladoga Lake from aboard a vessel with displacement 
of 120 tons (Fig. 8). The scope of tests included the measure- 


ments of the absolute value of gravity and DOV on moored 
vessel and during the vessel positioning at some points of the 
water area. 

Before installing the absolute gravimeter and the zenith 
telescope on board, gravity and DOV were measured at the 
pier on a massive concrete base. These measurements were 
further used as reference. To account for the difference in 
altitudes between the location of the gravity measurements 
on a concrete base and the place where the absolute gravime- 
ter was installed on the vessel, we performed leveling. The 
difference in the altitudes was 2.1 m. That was taken into 
account in the further data processing. 

The sea tests started from measuring the absolute value 
of gravity on the vessel moored at the pier. Within 1.5 h, 
a required data set of 100 reliable measurements taken at 
the inertial disturbances less than 3 mGal was formed. The 
difference between the gravity measurements at the pier and 
aboard the vessel was less than 0.1 mGal (Fig. 9). This 
result confirmed the possibility of offshore surveys without 
referencing to onshore gravity stations. 

The second stage of the sea tests was verification of 
gravity measurement accuracy when positioning the vessel 
at a point of the Ladoga Lake water area. The measure- 
ments were taken for 2 h at the sea state 2. The inertial 
accelerations reached 1.5 Gal, and the vessel pitching/rolling 
was up to 5°. In spite of significant accelerations (Fig. 10), 
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mGal 


Fig. 10 Gravity measurements at water area point 


25 reliable measurements were obtained. The difference of 
gravity values relative to the pier, calculated by the readings 
of the absolute and relative gravimeters, was +4.20 and 
+5.22 mGal, respectively. 

The obtained results confirmed the possibility of measur- 
ing the absolute value of gravity in real sea conditions at 
the sea state up to 2 points, using a system comprising the 
absolute and relative gravimeters. 

DOV measurements with the gyrostabilized zenith tele- 
scope were carried out in the water area of the Ladoga Lake, 
at two points located 100 m away from each other, on four 
nights. The results of the measurements were compared to 
the known values of DOV according to the global model 
EGM-2008, since they can be used as a reference in a 
region with small gravity anomalies and smooth topography 
(Fig. 11). 


min 


Wz (Chekan-AM) 
.  g(GABL-PM) 


Ni i 


min 


The standard deviations of DOV measurements were 
0.657 (Ł-component) and 0.85/77 (y-component). The differ- 
ences with EGM-2008 data were 0.967 (€-component) and 
0.36 (ņn-component). So the results of the field tests demon- 
strated the possibility of rapid high-precision measurement 
of the absolute values of DOV in marine conditions, using 
the astrogeodetic method. 


7 Conclusions 


The first results of experimental studies confirmed the possi- 
bility of high-precision measurement of the absolute values 
of gravity and deflection of the vertical at sea. The use of an 
integrated system will improve the quality of geodetic data 
in regions of the World Ocean with large gravity anomalies. 
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Fig. 11 The results of DOV measurements at water area point 
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Abstract 


In recent years, it was shown that the quality of strapdown airborne gravimetry using a 
navigation-grade strapdown inertial measurement unit (IMU) could be on par with “classi- 
cal” airborne gravimeters as the 2-axis stabilized LaCoste and Romberg S-type gravimeter. 
Basically, two processing approaches exist in strapdown gravimetry. Applying the indirect 
method (also referred to as “inertial navigation approach” or “one-step approach”), all 
observations — raw GNSS observations or position solutions, IMU specific force and angular 
rate measurements — are combined in a single Kalman Filter. Alternatively, applying the 
direct method (also referred to as “accelerometry approach” or “cascaded approach”), 
GNSS position solutions are numerically differentiated twice to get the vehicle’s kinematic 
acceleration, which is then directly removed from the IMU specific force measurement in 
order to obtain gravity. In the scope of this paper, test runs for the application of strapdown 
airborne and shipborne gravimetry are evaluated using an iMAR iNAV-RQH-1003 IMU. 
Results of the direct and the indirect methods are compared to each other. Additionally, a 
short introduction to the processing scheme of the Chekan-AM gravimeter data is given 
and differences between Chekan-AM and strapdown results of the shipborne campaigns are 
analysed. Using the same data set, the cross-over residuals suggest a similar accuracy of 
0.39 mGal for the Chekan-AM and 0.41 mGal for the adjusted strapdown results (direct 
method). 


Keywords 


Airborne gravimetry - Gravity disturbance - IMU - Shipborne gravimetry - Strapdown 


1 Introduction 


Airborne and shipborne gravimetry enables gravity determi- 
nation with a medium accuracy, compared to satellite and 
terrestrial gravimetry. Especially in remote areas, where ter- 
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Regions without satellite gravity data, it is the only practica- 
ble gravity determination approach. 

While shipborne gravimetry using gas-pressured, pendu- 
lum or spring gravimeters had already been employed in the 
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et al. 2005), first airborne gravity test flight results using 
a horizontally stabilized spring gravimeter have not been 
published before 1960 (Nettleton et al. 1960). In the 1990s, 
a major improvement in the obtained accuracy followed the 
advent of the Global Positioning System (GPS), followed 
by further Global Navigation Satellite Systems (GNSS). 
The application of inertial measurement units (IMUs) as 
gravimeters could be shown to be capable of generating 
results on an accuracy level similar to that of the “classical” 
approach, where horizontally stabilized spring gravimeters 
are used (Glennie et al. 2000; Becker et al. 2016). 

Main advantages of the IMU usage in the so-called 
“strapdown” approach are a lower instrument weight, less 
space and power requirements as well as lowered purchase 
and maintenance costs. Furthermore, the measurement prin- 
ciple enables vector (3-D) gravimetry and the system is 
less impaired by turbulences and flight manoeuvres. Strong 
sensor drifts affecting strapdown gravimetry can be reduced 
using calibration methods that can lower their impact signif- 
icantly. It was demonstrated that the application of a simple 
warm-up calibration of the IMU’s vertical accelerometer can 
eliminate most of the sensor drift (Becker et al. 2015b; 
Becker 2016). 

Even though the vehicle dynamics differ in airborne and 
shipborne gravimetry, processing strategies virtually agree. 

This paper gives an overview on strapdown gravimetry 
and recent campaigns with participation of the chair of 
Physical and Satellite Geodesy (PSGD) of the Technical 
University of Darmstadt. An introduction in the processing 
approaches in strapdown gravimetry is given in Sect. 2. To 
enable an evaluation of the presented strapdown result’s exte- 
rior accuracy, the data processing in “classical” kinematic 
gravimetry is briefly presented in Sect. 3 using the example 
of the Chekan-AM gravimeter of the German Research 
Centre for Geosciences (GFZ Potsdam). The airborne and 
shipborne campaigns evaluated in the scope of this paper 
are discussed in Sects. 4 (campaign introduction), 5 (strap- 
down processing results) and 6 (comparison to Chekan-AM 
results), followed by the concluding Sect. 7. 


2 Processing Approaches in Strapdown 
Gravimetry 


There exist two fundamentally different processing 
approaches in (kinematic) strapdown gravimetry: the indirect 
and the direct method. This paper shortly introduces some 
basics of both approaches; for a more detailed view, the 
reader is referred to Becker (2016) and Johann et al. (2019). 

The “indirect” method as it was introduced by Jekeli 
(2001) has also been called “inertial navigation approach” 
(Ayres-Sampaio et al. 2015), “traditional way” (Kwon and 
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Jekeli 2001) or “one-step approach” (Becker 2016). Here, 
all observations from an IMU and a GNSS receiver are 
integrated in a Kalman filter with a state vector including 
position, velocity, attitude, sensor biases and gravity. In 
the prediction step, the previous epoch’s state vector is 
propagated using the IMU observations (accelerations and 
angular velocities). If GNSS observations are available as 
pseudoranges and carrier phases (tightly coupled integration) 
or position and velocity solutions (loosely coupled integra- 
tion), the Kalman filter will form a weighted average between 
the predicted and the actual observations. The weighting 
depends on the accuracy ratio of the predicted states and 
GNSS observations as reflected by the stochastic model. 

“Accelerometry approach” (Ayres-Sampaio et al. 2015; 
Kwon and Jekeli 2001), or “cascaded approach” (Becker 
2016) have been alternative terms for the “direct” method. In 
contrast to the indirect method, where gravity determination 
is done in the position domain, in the direct method, vertical 
gravity is directly obtained in the acceleration domain. The 
gravity disturbance 

bg" = r” —Chf? + bgt, y" (1) 
in the navigation frame being the difference of gravity and 
normal gravity y” is directly obtained by subtracting the 
specific force C} f b from the kinematic acceleration #” 
(Wei and Schwarz 1998). All these terms are given in the 
navigation frame n whose origin is defined to be at the IMU’s 
centre of observations with its orthogonal axes pointing 
towards North, East and Down. The specific force f? as 
measured by the IMU accelerometers in the body frame b 
is transformed to the navigation frame by the multiplication 
with the rotation matrix C}. Alternatively, quaternions can 
be used to represent the attitude. The GNSS position solu- 
tions are numerically differentiated twice in order to obtain 
the kinematic acceleration. The vehicle’s attitude used as 
input for the rotation matrix is computed in a GNSS/IMU 
integration algorithm. The addition of the velocity-dependent 
Eötvös correction ô g” eliminates the impacts of the Coriolis 
acceleration and the centrifugal acceleration caused by the 
movement of the vehicle relative to the Earth. The obtained 
gravity disturbance needs to be low-pass filtered in order to 
eliminate high-frequency noise. 

On the one hand, the indirect method implements a well- 
known optimal estimation procedure, the Kalman filter, and 
is commonly deemed to be the more rigorous approach 
(Becker 2016). On the other hand, if the direct method 
is applied, the GNSS/IMU integration can be performed 
using available commercial navigation software (e.g. the 
NovAtel Waypoint InertialExplorer 8.60 used in the scope 
of this paper), significantly reducing the burden of software 
development. 
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In order to get the best out of the strapdown data, several 
corrections need to be applied. A detailed description can be 
found in Johann et al. (2019). 

e A bias and a linear drift are removed from the gravity 
disturbance results by anchoring the results to known ref- 
erence values at the airports/harbours combined with a lin- 
ear interpolation between these static periods. Usually, the 
reference values are obtained from local terrestrial gravity 
measurements of superior accuracy (about 0.05mGal = 
0.05 - 10->m/s? or better). 

e The distance between the IMU’s centre of observation and 
the GNSS antenna phase centre, the so-called “lever arm”, 
needs to be taken into account. 

e If no thermal stabilization housing is used, the gravity 
results of the iMAR iNAV-RQH-1003 IMU can be sig- 
nificantly improved with a thermal calibration. For the 
campaign analysis at hand, a warm-up calibration of the 
vertical accelerometer (Becker et al. 2015a) is applied. 

e In various strapdown campaigns, the authors have 
observed an heading-dependent shift of the gravity 
disturbance estimates. The effect is strongest when 
heading northwards or southwards (with opposite sign) 
and weakest when heading westwards or eastwards. 
In order to attenuate this systematic error, an empiric 
correction is applied: The cosine of the heading multiplied 
with a constant is added to the preliminary gravity 
disturbance results. The constant is selected empirically 
for the whole gravimetry campaign. The cause for this 
effect is unknown and is an object of research. Possible 
explanations might be an uncorrected instrumental error 
of the iMAR iNAV-RQH-1003 (Becker 2016) or a 
modelling/approximation error dependent on the current 
direction of travel. 

e Even though gravity disturbance is determinable during 
turns in strapdown gravimetry, the quality of the results in 
these track segments is lowered significantly. Therefore, 
to get accuracy indicators for the best available results and 
to be less dependent on the campaign trajectory design, 
only the data of approximately straight track segments 
(“lines”) is included in the quality analysis. 

In the scope of this paper, the precision of the results is 
examined using a cross-over analysis. For this, the survey 
plan includes intersecting lines. When the intersections are 
at nearly the same height and one assumes that there is no 
significant temporal gravity variation, the root mean square 
(RMS) of the cross-over residuals at these intersections is a 
measure for the precision of the results. Assuming an equal 
accuracy on both lines, the RMS error (RMSE), which is 
given by the RMS divided by V2, indicates the uncertainty of 
the obtained gravity disturbance values at the lines (Forsberg 
and Olesen 2010; Becker 2016). 

Since the gravity disturbance results have already been 
fixed to the reference values at the airports/harbours, the 


linear trend is easily removed from the results. If a sufficient 
number of cross-over points is available, the impact of further 
drifts during a measurement flight/cruise can be reduced 
by a cross-over adjustment (also called “levelling”). During 
the least-squares adjustment, a bias per line is estimated. 
If the number of cross-over points is low, the resulting 
precision value is susceptible to being too optimistic, which 
can be avoided by the use of correction factors (Becker 2016; 
Johann et al. 2019). 


3 Chekan-AM Data Processing at GFZ 


To enable a comparison with another measurement system, 
processing methods and some results using a Chekan-AM 
gravimeter are presented in this paper in addition to the 
strapdown analysis. 

Since 2011, the GFZ Potsdam has been involved in var- 
ious shipborne and airborne gravimetry campaigns with its 
mobile gravimeter Chekan-AM whose working principle is 
based on angle variation measurements of two quartz sensors 
that are positioned in a viscous liquid. The performance of 
this equipment has been verified in GFZ’s marine test cam- 
paign in Lake Miiritz, and other dedicated campaigns such 
as Lake Constance and the GEOHALO mission (Petrovic 
et al. 2015). The same equipment is further reliably used in 
the shipborne gravimetry campaigns in the Baltic Sea within 
the subactivity 2.1 of the “Finalising Surveys for the Baltic 
Motorways of the Sea” (FAMOS) project (Swedish Maritime 
Organisation 2019). 

The recovery of the measurements in terms of acceler- 
ations is based on the mathematical model and calibration 
constants provided by the manufacturer Electropribor, St. 
Petersburg. It is known that the survey vehicle and the mea- 
surement conditions have a strong impact on the gravimeter 
measurements. In principle, the processing scheme presented 
in Krasnov et al. (2011) is followed in GFZ’s data process- 
ing routine but the processing needs to be adopted based 
on the quality of the raw measurements and measurement 
conditions. Obviously, in mobile gravimeter measurements, 
the raw recordings need to be converted into acceleration 
units. Like in strapdown gravimetry, measurements need 
to be corrected for the Eötvös effect and, most impor- 
tantly, low-pass filtered to eliminate the high frequency 
noise components and retrieve meaningful results. Occa- 
sionally, different low pass filter lengths can be applied to 
different quality of datasets. Based on their experiences, 
GFZ found that, in shipborne gravimetry, a filter length of 
400 s provides satisfying results in terms of measurement 
accuracy w.r.t. global models and other campaign results 
computed in a cross-over analysis and is a good trade- 
off between the spatial resolution and final data accuracy 
(Lu et al. 2019). 


One of the most challenging tasks of relative gravimetry 
is the drift estimation of the gravimeter, which needs to be 
computed as good as possible and taken into account in the 
data processing. The gravimeter recordings at the harbours, 
where reference gravity values are available, are taken before 
and after completing relevant tracks and are used to compute 
the drift behaviour of the instrument sensors. After correcting 
the drift, since the measured values are not absolute gravity 
but gravity differences measured along the measurement 
track they need to be tied to absolute gravity values at 
reference stations that are ideally co-located to the harbour 
tie points. 

It is worth mentioning that the Chekan-AM measurements 
even after applying a low pass filter are not of very high 
quality during the turns of the ship. Therefore, these periods 
lasting about from few to tens of minutes as well as disturbed 
measurements due to harsh sea conditions or instrumental 
errors (e.g. unexpected drift behaviour, sensor aging) are 
manually detected and removed from the final delivered 
products. Good quality measurements (e.g. performed under 
optimum measurement conditions during dedicated gravime- 
try campaigns) may not require detailed investigations. 


4 Investigated Airborne and Shipborne 
Strapdown Gravimetry Campaigns 


Since 2013, PSGD co-organized or supported various strap- 

down gravimetry campaigns. An IMU of the type iMAR 

iNAV-RQH-1003 (Fig. 1, bottom right) (IMAR Navigation 

2012) was used as gravimeter. The IMU was connected to 

GNSS antennas for the purpose of time synchronization. The 

following campaign results are presented in this paper: 

— “MY2014”: In August 2014, 12 flights were undertaken 
above the South China Sea northwest of Borneo, 
Malaysia, using a medium-size aircraft, type Beechcraft 
King Air 350 (Fig.2a), with autopilot. For details, 


Fig. 1 Chekan-AM and iMAR iNAV-RQH-1003 with iTempStab- 
AddOn aboard DENEB 2018 
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including a scatter plot of the results and an elaborate 
evaluation, see (Johann et al. 2019). 

“ODW2017”: In May 2017, a small-scale survey was 
arranged with a motor glider of the type Grob G 109B 
(Fig. 2b). The flights were conducted at the Upper Rhine 
Graben above the Odenwald, a low mountain range with 
significant gravity variations, located in the southeast 
of Darmstadt, Germany. An objective of the flight was 
to test the performance of strapdown gravimetry under 
difficult conditions with a relatively small aircraft fly- 
ing in drape mode (following the topography’s altitude 
changes) without autopilot. Due to the small weight 
of the aircraft and the missing autopilot, the flight has 
been much less stable than MY2014. GNSS data acqui- 
sition was affected by vibrations, possibly also defor- 
mations, of the wing, where the antenna was placed 
on. 

“ODW2018”: In March 2018, a test flight was conducted 
above the same region with a Cessna 206 “Stationair 6” 
(Fig. 2c), a light aircraft of medium size compared to 
the aircrafts of MY2014 and ODW2017. Compared to 
ODW2017, the aircraft change allowed for a smoother 
flight, but sensor drifts caused by temperature changes 
had a greater impact in ODW2018 since temperatures 
before the start were very low, being near the freezing 
point. 

“BTS2017”: Among other campaigns, GFZ performed 
a dedicated survey in July 2017 on the German Fed- 
eral Maritime and Hydrographic Agency’s (BSH) survey, 
wreck-search and research vessel DENEB (Fig. 2d). The 
DENEB campaign measurement tracks were planned by 
the Federal Agency for Cartography and Geodesy (BKG) 
considering the gravity measurement densification needs 
for the determination of a high accuracy German quasi- 
geoid and the gravity measurements were collected by the 
GFZ’s Chekan-AM relative gravimeter. Additionally, in 
order to investigate the potential of strapdown shipborne 
gravimetry using navigation-grade IMUs, PSGD’s iMAR 
iNAV-RQH-1003 was installed aboard the DENEB side- 
by-side with the Chekan-AM. The geographical focus 
of the 2017 campaign was the Bay of Mecklenburg, 
and an area between the German island Riigen and the 
Danish islands Bornholm and Møn. Compared to airborne 
gravimetry, the Baltic Sea campaigns are challenging 
for strapdown gravimetry since some cruises (runs from 
harbour to harbour) are extraordinarily long (up to 60 h). 
“BTS2018”: In July/August 2018, the mission of the 
DENEB was continued with another campaign that 
focussed on the western (Bay of Kiel) and eastern areas 
(Bay of Pomerania including Polish territory) of the 
German part of the Baltic Sea. In comparison to the 
2017 campaign, the strapdown measurement system was 
upgraded by encasing the iMAR iNAV-RQH-1003 in a 
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Fig. 2 Vehicles: (a) Beechcraft King Air 350 (MY2014); (b) Grob G109B (ODW2017); (c) Cessna 206 “Stationair 6” (ODW2018); (d) German 
research vessel DENEB (BTS2017/18) 


Table 1 Strapdown gravimetry campaign overview (MY2014 results taken from Johann et al. (2019)) 


BTS2017 BTS2018 
Aircraft/vessel Cessna 206 “Stationair 6” DENEB 

Number of fighis/eruses P EEE ise E Ie 

Total Tine distance [km] 25 
Mean line velocity [m/s] 48 


Ellipsoidal flight height 1953/1016/4256 655/396/927 | 947/797/980 
(mean/min/max) [m] 


Number of cross-over points |0 œ ë a ë ë a ë oo 
Turbulence (RMS-g) [mm/s?] 48 
GNSS sampling frequency z [5 aoo o o a 
Low-pss er e R e E E e Eom 
Half-wavelength resolution [km] 


RMSE non-adjusted [mGal] direct | 1. 26 (1.3) 3. 10 (3.0) 3. 97 1.28 (1. C a 81 7 H 
(indirect) 

RMSE adjusted [mGal] direct 0.62 (0.68) 1.12 (0.90) 1.00 (0.81) | 0.86° | 0.55° 
(indirect) 


*BTS2017: Three cruises are remaining if one of the four cruises was rejected assuming an outlier 

>BTS2017: The heading-dependent correction has not been applied in the indirect method 

©The cross-over adjustment of the BTS campaigns was performed stint-wise (complete run from harbour to harbour), while the other campaigns 
were adjusted line-wise. Usually, lower RMSE values are obtained via line-wise adjustment 


housing, the iMAR iTempStab-AddOn, which stabilises 


the IMU’s temperature using two Peltier elements. 60 
Dual frequency GNSS receivers tracking GPS and _ 
GLONASS were installed in all campaigns. Positions were = aa 
calculated in the Precise Point Positioning (PPP) mode with E ea 
the NovAtel Waypoint InertialExplorer 8.60 using precise w 
satellite orbit and clock products by the Center for Orbit 0 
Determination in Europe (CODE). The campaign statistics 
are summarised in Table 1. The turbulence metric and the -20 
obtained precision values (RMSE) will be discussed in 8.5 8.6 8.7 8.8 8.9 9 
Sect. 5. Longitude [°] 


Fig. 3 Vertical gravity disturbance [mGal] ODW2017 (Map data: 
Google, GeoBasis-DE/BKG, Europa Technologies) 


5 Strapdown Results 


The obtained vertical gravity disturbance without crossover campaign, equivalent figures including preliminary results of 
adjustment for both Odenwald and both Baltic Sea cam- the horizontal components can be found in (Johann et al. 
paigns are illustrated in Figs. 3, 4, 5, 6. For the Malaysia 2019). RMSE obtained in cross-over analyses with and 
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Fig. 4 Vertical gravity disturbance [mGal] ODW2018 (Map data: 


Google, GeoBasis-DE/BKG, Europa Technologies) 
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Fig. 5 Vertical gravity disturbance [mGal] BTS2017 (Map data: 
Google, GeoBasis-DE/BKG, Landsat/Copernicus) 
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Fig. 6 Vertical gravity disturbance [mGal] BTS2018 (Map data: 
Google, GeoBasis-DE/BKG, Landsat/Copernicus) 


without cross-over adjustment for all four campaigns are 
noted in Table 1, where available results of the indirect 
method of strapdown gravimetry are given in parentheses. 
For the Malaysia 2014 campaign, the results of the direct 
and the indirect method are on a par with an RMSE of 
about 1.3 mGal without and 0.65 mGal with line-wise cross- 
over adjustment. For the Odenwald 2017 campaign (Fig. 3), 
results are considerably worse with an RMSE of 3.1 mGal 
without and 1.1 mGal with adjustment. The indirect method 
performs slightly better (3.0 and 0.9 mGal). Without adjust- 
ment, the results for Odenwald 2018 (Fig. 4) are even worse 
(4.0 mGal), maybe caused by temperature calibration prob- 
lems due to the low environment temperature. In contrast, 
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after adjustment, the RMSE is at the precision level of the 
Malaysia campaign (0.64 mGal). 

For the Baltic Sea 2017 campaign (Fig. 5), without adjust- 
ment, an RMSE of 1.28 mGal is obtained applying the 
direct method and 1.18 mGal applying the indirect method. 
It should be noted that the empirical heading-dependent 
correction was not used in the indirect method. Hence, 
further improvements might be possible. After a cruise- 
wise cross-over adjustment, RMSE improved to 1.00 mGal 
(direct method) and 0.81 mGal (indirect method). A line- 
wise adjustment was not possible, since isolated nets would 
occur causing a datum defect. One of the four cruises has 
been identified as an outlier. The rejection resulted in an 
improved accuracy of 0.81 mGal without adjustment (direct 
method). A cruise-wise adjustment did not enhance the 
results after the outlier rejection. In the Baltic Sea 2018 
campaign (Fig. 6) with the installed temperature stabilisation 
housing, significantly better results with an RMSE of 0.71 
mGal before and 0.55 mGal after the cruise-wise adjustment 
have been noted. Again, a line-wise adjustment was not 
applicable. To test the repeatability of the results, both Baltic 
Sea campaigns have also been examined together (without 
the outlier cruise in 2017) leading to RMSE of 0.87 mGal 
without and 0.66 mGal after adjustment, medium accuracies 
compared to the campaign’s single results. 

The accuracy of airborne gravimetry results is suspected 
to depend on the flight turbulence level. A turbulence metric 
that directly indicates the reaction of the aircraft turbulences 
is the RMS-g. It is defined as the moving standard deviation 
of the vertical acceleration during the flight (Becker 2016). 
Other turbulence metrics like the Eddy Dissipation rate, 
which rate the atmospheric turbulence state, eliminate some 
aircraft-dependent effects acting on the IMU. In the scope of 
this paper, the RMS-g is computed applying a 50 s moving 
window on | Hz GNSS vertical acceleration data. For the 
presented campaigns, a trend between cross-over RMSE and 
RMS-g might be indicated by the results given in Table 1 
if the non-adjusted RMSE of ODW2018 was declared as an 
outlier caused by temperature calibration issues. Note that 
accuracy differences depend on multiple campaign parame- 
ters like the local gravity field, the low-pass filter length and 
GNSS data quality. 


6 Comparison to Chekan-AM Results 


Using the cross-over analysis algorithm implemented in the 
direct method of strapdown gravimetry, an RMSE of 0.33 
mGal and 0.39 mGal resulted for the Chekan-AM data of the 
Baltic Sea campaigns 2017 and 2018, respectively. To enable 
a comparison, line segments affected by harsh sea conditions 
have also been removed from the strapdown results resulting 
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Fig. 7 Comparison of estimated vertical gravity disturbances for BTS2018, cruise 2: (a) blue: strapdown (direct), green: Chekan-AM; (b) 


Difference (RMS: 0.9 mGal) 


in a slightly improved RMSE for BTS2017 (improvement 
of about 0.1 mGal compared to Table 1). For BTS2018, the 
strapdown RMSE decreased to 0.52 mGal before and 0.41 
mGal after cruise-wise adjustment. Regarding the absolute 
differences to the non-adjusted strapdown results (direct 
method), the RMS are 1.24 mGal in 2017 (without the outlier 
cruise) and 0.97 mGal in 2018. This indicates remaining 
systematic differences, possibly due to different sensor drift 
behavior. The comparison of the obtained vertical gravity 
disturbances of cruise 2, BTS2018, are shown in Fig. 7. 


7 Conclusions 


After applying several corrections, the precision of the pre- 
sented strapdown gravimetry campaigns is between 0.7 and 
4.0 mGal without cross-over adjustment. The best RMSE for 
an airborne campaign (1.26 mGal) was reached at Malaysia 
2014, which was conducted with a medium-size aircraft with 
autopilot under relatively stable flight conditions. 

After a line-wise adjustment, the airborne RMSE obtained 
at Malaysia 2014 and Odenwald 2018 were at the same 
level (about 0.6 mGal). The corresponding precision of the 
Odenwald 2017 campaign is significantly worse (1.1 mGal), 
probably suffering from GNSS acquisition problems. 

The successful combination of both Baltic Sea cam- 
paigns indicates a good repeatability of the direct method 
of strapdown gravimetry. In the implementations created 


at PSGD, the accuracy of this method is on a par with 
indirectly obtained results for most evaluated campaigns. 
Many additional airborne and shipborne surveys processed 
with both methods would be required to enable statistically 
firm statements about the accuracy ratio of both approaches 
under various conditions. 

The first test of PSGD’s new thermal housing at the Baltic 
Sea 2018 campaign indicates that its installation improves 
the long-term stability significantly. Comparing line seg- 
ments without harsh sea conditions, the precision of the 
cruise-wise adjusted strapdown gravity disturbance results 
of 0.41 mGal is competitive with the Chekan-AM results 
(0.39 mGal). Remaining systematic differences between the 
strapdown and Chekan results are a subject of research. 

Advantages of the strapdown approach like the possibility 
of vector gravimetry let this approach be an auspicious alter- 
native to the classical approach using horizontally stabilized 
spring gravimeters. 
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Abstract 


In this article, we propose a technique for processing airborne gravimetry measurements 
at repeat drape lines. To model gravity along a line, we use cubic B-splines. The unknown 
parameters of the model (the spline coefficients) are to be estimated from the measurements. 
In order to take into account dependance of gravity on position, the same set of the spline 
coefficients is assumed for gravity modeling at each repeat line. Gravity estimation is 
performed with the Kalman filter. The spline coefficients are included in the state vector of 
the filter as unknown constants. The developed estimation algorithm was tested using mea- 
surements of the GT-2A gravimeter collected at repeat drape lines. The gravity estimates 
obtained with the use of B-splines are more accurate compared with those obtained with the 
standard technique, in which gravity is modeled as a stochastic process in the time domain. 
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Airborne gravimetry - Gravity disturbance - GT-2A gravimeter - Kalman filter - Repeat 


drape lines - Splines 


1 Introduction 


Airborne scalar gravimetry is widely used to obtain high 
accuracy information about the Earth gravity in hard-to- 
reach regions (shelves, mountains, polar areas, etc.). A 
typical platformed airborne gravimetry system consists of 
a gravity sensor (a single-axis accelerometer) mounted on a 
gyro-stabilized platform and two (or more) global navigation 
satellite system (GNSS) receivers (one on board the aircraft, 
others on the ground). Postprocessing of airborne gravimetry 
measurements includes three main stages such as determina- 
tion of GNSS positioning and velocity solutions, estimation 
of the platform tilt angles via the inertial navigation system 
(INS) and GNSS integration, and estimation of gravity using 
the Kalman filtering technique or low-pass filtering. Both 
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filtering techniques are based on representing gravity in the 
time domain (e.g., (Lu et al. 2017) or modeling gravity as 
a stochastic process (Bolotin and Golovan 2013; Kwon and 
Jekeli 2002; Stepanov et al. 2015)). 

In this research, we focus on the last stage of 
postprocessing (i.e., gravity estimation) given airborne 
measurements at repeat drape lines. Using a stable platform 
gravimetry system (e.g., the GT-2A (Berzhitzky et al. 2002)) 
in a drape flight is challenging, since such a system is more 
suitable for benign dynamics (Studinger et al. 2008). For the 
GT-2A, underestimation of the scale factor of the gravity 
sensor can be observed during postprocessing of drape lines. 
For this reason, specific calibration lines are flown, at which 
the scale factor is reliably estimated. However, accuracy of 
gravity estimation at repeat drape lines is still lower than in 
the case of a flight at a constant height. To our opinion, the 
standard approach to gravity estimation based on modeling 
gravity in the time domain is not well suitable for processing 
repeat lines, as gravity modeled along a line, strictly 
speaking, differs from gravity modeled along a repeat line. 
Actual gravity is, of course, the same at repeat lines assuming 
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that they were flown at the same altitude. Therefore, in the 
case of a repeat line flight, it is preferable to model gravity 
in the spatial domain rather than in the time domain. 

In this work, we propose an approach based on spatial 
modeling of gravity to process repeat drape lines flown 
over the same ground track. We parameterize the gravity 
disturbance along a line using B-splines of order three. To 
take into account dependance of gravity on position, we use 
the same set of unknown parameters of the model (the spline 
coefficients) for gravity modeling at each repeat line. The 
model is easily described in the state space assuming that 
the spline coefficients are constant over the time during the 
aircraft’s flight. The standard stochastic estimation problem 
is posed, and solved via Kalman filtering. The state vector for 
the filter includes the spline coefficients and the parameters 
of the gravimeter systematic errors (the calibration parame- 
ters). 

The developed approach was applied to processing air- 
borne measurements collected with the GT-2A system during 
a repeat drape line flight in South Africa in 2009. The 
obtained gravity estimates are compared with those provided 
by the standard approach based on modeling gravity in the 
time domain, and with the gravity values derived from the 
global model EGM2008. An increase in accuracy in gravity 
estimation using the new approach is shown. 


2 Gravity Estimation Methodology 


Basic Equation and Gravimeter Measurements The grav- 
ity sensor of the GT-2A is mounted on a gyrostabilized 
platform, which makes the instrument sensitivity axis to be 
vertical. The gravity disturbance ôg is determined from the 
fundamental equation of airborne scalar gravimetry, which is 
written as: 

V3 = —ôg — go + gen + f (1) 
where V3 is the vertical velocity of the gravimeter proof mass 
in the ENU coordinate frame, go is the normal gravity at the 
flight height, g¢;, is the Eötvös correction term, f3 is the 
vertical projection of the specific force. 


Raw measurements of the gravimeter fy can be described 
by the equation 


A +k) fi + tft = fa — ô fhor + Ofritr + 6 farift + 4f- 


Here k3 is the scale factor error, t is the time delay of 
gravimeter measurements (due to physical delay and aver- 
aging of gravity sensor measurements performed in the 
internal software of the gravimeter), fhor is the horizontal 
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acceleration influence on fy due to imperfect placement 
of the gravimeter on the platform, i.e., more explicitly, 
Sfhor = ki fl + ko ff, where ff, f3 are the platform 
accelerometer measurements, kı, k2 are unknown factors that 
can differ from one flight to another. Further, 6 friis is the 
horizontal acceleration due to the platform tilt error (it is 
determined from INS-GNSS integration), 6 farift is the sum 
of the gravity sensor bias and drift. Assuming linear drift, 
the term 6 farift can be computed given pre- and post-flight 
measurements and terrestrial gravity value. The term qp is 
the noise. 

The frequency of the GT-2A raw measurements is 
18.75 Hz, and the time delay t is about 2 s. 

As the terms 6 frin and 6 farift can be determined and 
compensated for using standard procedures, we will omit 
them further on. 


State-Space Modeling Rewrite Eq.(1) replacing unknown 
variables on the right-hand side of the equation by measure- 
ments and omitting the gravity disturbance term: 


Vs = -g0 + 25 + fh (2) 


where go and gz, are known from the GNSS position and 
velocity solutions (assuming that the lever arm effect is 
compensated for). 


Denote by AV; the vertical velocity error defined as 
(Bolotin and Golovan 2013) 

AV = V3 = V3 =T fy. (3) 

Then we obtain the error dynamics equation of airborne 

scalar gravimetry by subtracting Eq.(2) from Eq.(1) and 
using the above gravimeter measurement model 


AV; = —ôg + kifl thf +k fi -qr (4) 


where kı, k2 are the residual platform misalignment errors 
in the sum with errors of placement of the gravity sensor on 
the platform, and k3 is the scale factor error of the gravity 
sensor. 

Representing the GNSS velocity solution as V; 
V3 + vei, where eye; is the GNSS velocity random error, and 
introducing the velocity observation as y = yor SS _ V3, we 
obtain the observation equation 


GNSS _ 
VENSS = 


y=AV34+ Tt ff + ever. (5) 
Equations (4) and (5) contain six unknown variables: 


AV3, bg, T, ky, ko, k3. 


The variables q f and eve; in Eqs. (4) and (5) are assumed to 
be zero-mean white noises with known variances. We adopt 
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the following models for the gravimeter systematic error 
parameters (calibration parameters): 

T=4r, ki = dq. ko=de, k3 = qk (6) 
where qr, 4k,»>9ko, 4k; are zero-mean white noises with 
known (sufficiently small) variances. 

By introducing an equation for dg (a gravity model), we 
obtain the system of equations in the state-space form. The 
gravity estimation algorithm includes the following steps: 

1. Determining Vj from Eq. (2). 

2. Estimation of gravity and the gravimeter calibration 
parameters via Kalman filtering and smoothing given 
the state-space system Eqs.(4)-(6) and the equation 


describing a gravity model. 


3 Modeling Gravity Along Repeat Lines 
Using B-Splines 


To arrive at estimation algorithm, we introduce a spatial 
gravity model assuming a repeat line flight. We use B- 
splines of order three for gravity modeling, which form 
a basis in the space of cubic splines in R! and each 
function of which has finite support (Fig. 1). Recall that 
the B-spline, by definition, is a fourfold convolution of 
a rectangular function with itself (De Boor 2001). We 
assume that the repeat lines are flown along the same 
ground track at the same altitude and have equal length 
L. Let x = x(t) € R! denote the distance travelled by 
the aircraft (more precisely, the gravity sensor proof mass) 
along the first line. Then x varies from 0 to L. For simplicity, 
assume that the subsequent repeat lines coincide with the 
first line. Thus, x varies from 0 to L or from L to O at 
the repeat lines, depending on the direction of aircraft’s 
flight. 

We represent the gravity disturbance at a point of a line as 
follows 


N-1 
5g(x) = D> ci Bi(x) (7) 
i=0 


Fig. 1 Cubic B-splines 


where B;(x) is the cubic B-spline, c; is the i-th spline 
coefficient to be estimated, N is the total number of the 
spline coefficients. Note that the same spline coefficients 
ci, i = 0,...,N — 1, are used for describing gravity at 
each repeat line. Note also that the B-spline B;(x) in fact 
is the function of (x — x) /Ax, where x T ica are the 
knots located at the first line, Ax is the length of the knot 
step (measured along the line). The value of Ax is to be 
chosen. 

Assuming that the spline coefficients are constant over 
the time during the flight, we can easily transform the 
gravity model Eq. (7) into the state-space form by writing the 
obvious equations 


(8) 


Therefore, we obtain the complete state-space system 
described by Eqs. (4)—(6) and (8) given the gravity model (7). 
Using Kalman filtering, the optimal (in the minimum mean 
squared error sense) estimates of the spline coefficients, 
Ci, i = 0,...,N — 1, the velocity error AV3 and the 
calibration parameters, T, kı, k2, k3, are obtained. The total 
dimension of the state vector is N + 5. At the initial 
iteration of the filter, the a priori variance of the state vector, 
including the variance of the spline coefficients, should be 
chosen. 


Transfer Function of the Optimal Filter The Fourier trans- 
form of a cubic B-spline can be written as (De Boor 2001): 


7 ef eT —] 
Bo(@) = (<>) 
joT 


where T = Ax/V is the time step of the knots, V is 
the average speed of the aircraft at the lines, j = V—I. 
By varying the time step of the knots T, one can obtain 
B-splines with different frequency bandwidths. Since then, 
it is desirable to determine how the spatial resolution of 
the gravity estimate provided by the Kalman filter depends 
onT. 


4 


For this, let us derive the expression for the transfer 
function of the filter. Consider the state-space equations 
Eqs. (4)-(8) in the discrete-time form denoting by Af the 
measurement time step. Neglect for simplicity the gravimeter 
calibration parameters and the gravity sensor noise qy in 
the equations. Let the measurement time step Ar and the 
spline knot step T satisfy the relationship T = mAt for 
some integer m > 2 and assume no a priori information 
on the variances of the spline coefficients (i.e., the variance 
of each coefficient tends to infinity). As above, assume that 
the GNSS velocity error that appears in Eq. (5) is a white 
noise. Then it can be shown that the optimal filter transfer 
function which maps the acceleration observations to the 
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Fig. 2 Transfer functions of the optimal filter in the spline-based 
approach and the standard approach based on the Gauss-Markov model 
of order 2 for gravity, and the Butterworth filter (two-pass) of order 2 


gravity disturbance estimate can be written as (Unser et al. 
1993) 


m—1 -1 
H(z) = LAG (+ Bice) 
k=0 


where Bo(z) is the Z-transform of the B-spline of order 
4. 

The obtained transfer function is close to the Butterworth 
filter (two-pass) of order 2 (see Fig. 2). Figure 2 also shows 
the transfer function of the optimal filter in the standard 
approach to gravity estimation, which is based on stochastic 
modeling of gravity (e.g., using a Gauss-Markov process of 
order 2 with infinite correlation time, that is, a process whose 
second time derivative is a white noise). In this case, the 
transfer function of the filter is close to the Butterworth filter 
(two-pass) of order 4. 

Thus, with the use of the above expression for the transfer 
function of the optimal filter (for the spline-based approach), 
the relationship between the filter cutoff wavelength A, and 
the length of the spline knot step Ax (or the time step of the 
knots T) can be derived as 

Ag © 34x, 


Ax = VT. (9) 


4 Results 


Survey Flight Overview The developed approach based 
on spatial gravity modeling using B-splines was applied 
to airborne gravimetry data collected within a flight 
over the Vrederfort Dome impact crater (120km south- 
west of Johannesburg, South Africa) on March 21, 2009. 
The GT-2A gravimetry system used in this survey was 
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Fig. 3 Survey area and the ground track of the repeat line flight 


installed in a BN-2T Islander aircraft operated by Xcalibur 
Airborne Geophysics (Olson 2010). Two dual-frequency 
GPS receivers operating at 10 Hz were used during the flight 
(the rover on board the aircraft and the base placed on the 
ground). 

The flight consisted of six repeat drape lines (with the 
ID labels RO1—RO6) and two calibration lines (with the ID 
labels RO7—RO8). The latter were flown in order to obtain 
more accurate estimate of the gravimeter scale factor error 
k3. All the lines were flown in the north-south and south- 
north direction over the same ground track (Fig. 3). The drape 
lines were flown at the heights from 1,400 to 1,600m above 
the WGS-84 ellipsoid and accurately followed the terrain 
(Fig.4). The calibration lines were flown at the height of 
3,000 m (Fig. 5). 

The length of each line is 85km. Average flight speed 
is 57 m/s with standard deviation (STD) 1.3 m/s. The flight 
duration is 3h 50 min. 
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Postprocessing Strategy Postprocessing was performed Table 1 Estimates of the gravimeter calibration parameters: time 


using the GTNAV/GTGRAV suite of programs developed 
by the Navigation and Control Laboratory of the MSU 
(Bolotin and Golovan 2013) and included the following 
main steps: 
1. Computation of the GPS position and velocity solutions 
in the carrier phase differential mode. 
2. Estimation of the platform tilt angles via INS-GNSS 
integration. 
3. Computation of the normal gravity and Eötvös correction, 
and solving Eq. (2). 
4. Estimation of the gravimeter calibration parameters and 
the gravity disturbance along the flight path. 
At the last step, gravity was modeled in the time domain as 
a Gauss-Markov process of order 2 with infinite correlation 
time. 


Accuracy of Estimating the Scale Factor Error k3 The 
maximum vertical acceleration at the repeat drape lines 
(with the Eötvös effect and the normal gravity subtracted) 
equals 0.1 m/s? (Fig.6). Hence, for gravity estimation 
with the accuracy of 0.5mGal the scale factor error k3 
should be estimated with the accuracy of 0.00005 or 
better. 


delay t (s), misalignment errors kı, kọ (arcmin), and scale factor error 
k3 


T 
= 
= 


Table 2 Standard deviations of the estimate errors for t (s), kı, ko 
(arcmin), and k3 


Drape lines 
Calibration lines 


Line 
Drape lines 
Calibration lines 


0.00001 
0.021 [0.00001 


Table 1 shows the calibration parameter estimates 
obtained with the developed spline-based approach from 
processing the drape lines and the calibration lines. The 
corresponding standard deviations of the estimate errors 
provided by the Kalman filter are shown in Table 2. Equal 
estimation accuracy of 0.00001 for k3 is achieved at the 
drape lines and at the calibration lines. This accuracy is 
sufficient for the gravity estimation error due to k3 to be at 
the level of 0.5 mGal. Moreover, since the two estimates of 
the scale factor error k3 are close enough to each other (the 
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Fig. 6 Vertical acceleration smoothed with the cutoff frequency of 1/100 Hz 
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Fig. 7 The along-track gravity disturbance estimates obtained with the two approaches for the cutoff wavelength of (a) 100s and (b) 70s, and 
gravity from EGM2008; (c) GPS height (averaged over six repeat drape lines) 


difference equals 0.00005), it is sufficient to process only 
the repeat drape lines to obtain a reliable estimate of k3, and 
there is no necessity to use the calibration lines. It should 
be noted here that the standard approach, which is used in 
the GTGRAV software and based on modeling gravity in the 
time domain, provides a reliable estimate of k3 only from 
processing the calibration lines. 


5 Discussion 


The gravity disturbance was estimated with the new approach 
based on spatial modeling of gravity using B-splines and 
with the standard approach based on modeling gravity in 


the time domain. Figure 7 shows the along-track gravity 
disturbance estimates obtained for the filter cutoff wave- 
length (inverse cutoff frequency) of 100 and 70s, which 
are equivalent to the spatial resolution Ag/2 of 3 and 2km, 
respectively. From Eq. (9), the corresponding values of the 
length of the spline knot step Ax are 2 and 1.4km. These 
values resulted in N = 46 and N = 64 spline coefficients in 
the gravity model Eq. (7), respectively. 

Figure 8 shows how the estimates of the spline coefficients 
evolve during processing of the repeat drape lines. The 
estimate changes occur due to the measurement update step 
of the Kalman filter. Since B-splines have finite support in 
the time domain, the estimate changes are local. As the spline 
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Fig. 8 Evolution of the estimates 
of two spline coefficients over the 
time at six repeat drape lines. The 
vertical coloured stripes indicate 
the time intervals of the aircraft 
turns 


coefficients depend on position, each new repeat line slightly 
refines the coefficient estimates. 

Mean value of the difference of the gravity estimates 
obtained with the two approaches (100s cutoff wavelength) 
is 0.05 mGal, the STD is 1.13 mGal. It can be seen from 
Fig. 7 that the gravity estimates provided by the new algo- 
rithm are more detailed and better correspond to the height 
profile in the interval from 60 to 70 km. 

In the absence of ground gravity data, we compare the 
gravity estimates with the global gravity model EGM2008, 
which is sufficiently accurate in the area of interest, as 
dense gravity datasets were available at the model construc- 
tion stage. Table 3 summarizes the statistics for comparing 
the gravity estimates and the gravity values derived from 
EGM2008. It can be seen that both approaches (spatial mod- 
eling and modeling in the time domain) show the same result 
for the cutoff wavelength of 180s. For 320s (equivalent to 
the nominal spatial resolution of EGM2008, which equals 
9 km), the new approach showed greater discrepancy than the 
standard one. The reason is likely in short-wavelength errors 
present in the gravity estimate provided by the new approach, 
as the optimal filter has less steep slope than in the case of the 
standard approach (see Fig. 2). 

At last, the gravity estimates were computed by process- 
ing a different number of repeat drape lines (from 1 to 
6) to analyze how each new repeat line affects the grav- 
ity estimate. The estimates were obtained with the two 
approaches for the cutoff wavelength of 100s and compared 
to EGM2008. The results are shown in Fig.9. The depen- 
dance on the number of the processed lines is different for 
the two approaches. For the new approach, the discrepancy 
with EGM2008 grows with the number of processed lines, 
which can be explained by the more and more detailed 
along-track estimate of gravity. For the standard approach, 
on the contrary, the along-track gravity estimate obtained 
by averaging over the lines becomes less detailed with 
the number of processed lines, and tends to gravity from 
EGM2008. 


37 38 39 4 414 42 43 44 45 46 
Time (s) x 10° 


Table 3 Standard deviations of the difference of the along-track 
gravity estimates obtained with the two approaches (spline-based 
modeling and Gauss-Markov modeling) and gravity values derived 
from EGM2008 (mGal) 


Cutoff wavelength | Spatial modeling Modeling in 
(spatial res.) of gravity (splines) the time domain 
70s (2km) | 4.10 | 5.00 
100s (3 km) 4.24 4.99 
180s (5 km) 4.05 | 4.03 
320s (9 km) | 3.24 | 2.79 
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| |] 
45 
© 
oO 
= 


Gauss-Markov 
—— Splines 


1 2 3 4 5 6 
Number of processed repeated lines 


Fig. 9 Standard deviations of the difference of the gravity estimates 
obtained for the cutoff wavelength of 100s and gravity from EGM200 


6 Conclusions 


In this research, an approach to gravity estimation at repeat 
drape lines was developed using B-splines for spatial mod- 
eling of gravity. This approach was compared with the 
standard one based on modeling gravity in the time domain 
by processing the GT-2A gravimeter data. The spatial gravity 
modeling showed an increase in accuracy of estimating the 
gravity sensor scale factor significantly. The achieved accu- 
racy was sufficient for gravity estimation along the repeat 
drape lines with the accuracy of 0.5mGal. Moreover, the 
calibration lines, which are necessary for accurate estimation 
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of the scale factor when using the standard approach, are not 
required in the case of the developed approach. The reliable 
estimate of the scale factor was obtained from processing the 
repeat drape lines. 

The gravity estimate provided by the new approach is 
in good agreement with EGM2008 (for the filter cutoff 
wavelength close to the minimum wavelength represented 
by the global model). For the shorter cutoff wavelengths, the 
new approach provides more detailed gravity estimates than 
the standard one. 
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Abstract 


Marine gravimetric studies were carry out in the waters of the Sea of Japan and the 
Tatar Strait. Maps of the gravitational field of the earth were completed, structural density 
modeling was executed, the depth of the Mohorovičić discontinuity was determined. 
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1 Introduction 


The Sea of Japan, washing the southeastern coast of Russia, 
cuts almost across the cut the southern tip of the Sikhote- 
Alin and Laoelin-Grodekov fold systems. Here, at a short dis- 
tance, a radical restructuring of the earth’s crust takes place: 
the transition from a mature continent to a young oceanic 
crust with the disappearance or significant processing of 
the upper sialic shell rich in ore and non-metallic minerals. 
According to modern concepts, this, like the formation of 
the Sea of Japan as a whole, is the result of Meso-Cenozoic 
destruction of the margin of the Asian continent and rift- 
ing, which was replaced by spreading with the formation 
of a young oceanic crust in the eastern part of the deep- 
sea Japanese Sea basin. As a result, a region was formed 
within which two radically different types of the earth’s crust 
closely “coexist”. 

The deep structure of the earth’s crust of this region is still 
an actual object of study. The decisive role in this is assigned 
to geophysical methods, in particular, seismic sounding and 
gravimetry. The capabilities of gravimetry to study the prob- 
lem under consideration are based on the fact that one of 
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Mohorovičić discontinuity - Sea of Japan - 


the parameters distinguishing the continental crust from the 
oceanic crust is their significantly different thickness (the 
depth of the Mohorovičić discontinuity or Moho (M)). This 
parameter with acceptable reliability is determined using 
gravimetry in combination with the reference data of seismic 
sensing, which makes it possible to determine the relief of 
the base of the earth’s crust in a given area. On the other 
hand, with a “fixed” structural framework, density modeling 
of the geological environment is possible, which makes it 
possible to study changes in the material composition of 
the crust masses. The combination of both situations makes 
it possible to realize structural-density (structural-material) 
modeling and, thus, to study the deep structure of the earth’s 
crust at the junction of its heterogeneous types. 


2 Marine Gravimetric Studies 


The first measurements of gravity in the Sea of Japan were 
made by Soviet researchers in 1937. Professor of Moscow 
State University L.V. Sorokin in a submarine made measure- 
ments of gravity at 74 points (Stroev et al. 2007). 

Until the mid-1980s of the last century, the main con- 
tribution to the study of the anomalous gravitational field 
of the Sea of Japan was made by the P.K. Sternberg’s 
State Astronomical Institute, who is part of Moscow State 
University (GAISH Moscow State University), in particular, 
his staff V.L. Panteleev, A.D. Gainanov, P.A. Stroyev and 
others. Unfortunately, the measurements were performed on 
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disconnected polygons and profiles. Profile measurements 
were taken over the Pervenets Rise, over the southern part 
of the Bogorov Ridge, detailed areal observations within the 
Yamato Rise, at three polygons: near the Korea peninsula, 
and near the mainland slope of northern Primorye (Belousov 
et al. 1973; Stroev and Kovylin 1973). 

At the same time, the entire eastern half of the Sea of 
Japan was covered by gravimetric surveys by the Hydro- 
graphic Service of the Ministry of Defense of Japan and 
Japanese universities (Tomoda 1973; Tomoda and Fujimoto 
1981). 

Since the mid-1980s, the V.I. IV’ichev Pacific Oceanologi- 
cal Institute Far Eastern Branch Russian Academy of Science 
(POI FEB RAS) has concentrated geophysical work in the 
northwestern part of the Sea of Japan (Fig. 1), within the 
economic zone of the Russia (formerly USSR). In addition 
to onboard gravimetry, seismic profiling and magnetometry 


Fig. 1 Scheme of the 
gravimetric study of the Sea of 
Japan. OP54—54 cruise R/V 
“Academician Oparin”, LV81-81 
cruise R/V “Academician M. 
Lavrentiev”, OP55—55 cruise R/V 
Academician Oparin, LV85-85 
cruise R/V “Academic M. 
Lavrentiev”, GA — cruises made 
on vessels of the Far Eastern 
Branch of the Russian Academy 
of Sciences in the period from 
1982 to 2010 


Vladivostok 


Me 


were included in the set of these works (Prokudin et al. 
2018). 

The first measurements of gravity at this stage (1987— 
1989), simultaneously with continuous seismic profiling 
(NSP) and hydro-magnetic surveying, were performed using 
on-board gravimeters GAG-ZhZ and GMN-K with analog 
recording of recordings on recorders. At that time, the 
navigational reference was very weak, positioning in the 
open sea was carried out no more than once in 4 h, in coastal 
areas, data from coastal radio navigation stations were used. 
The coordinates of the observation points were calculated by 
interpolation and navigational calculations, location errors 
introduced serious errors in the final accuracy of the survey. 
The measurement error was + (2.4-3.0 mGal). 

In 1990 (the 6th cruise of the R/V “Professor Gagarin- 
sky”), survey was carried out on the border of the economic 
zones of Russia and North Korea, as well as on the Yamato 
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Pacific 
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Rise. The survey was carried out by six GMN-K instruments 
with analog registration (error + 2.4 mGal). The direction of 
the profiles is meridional on the Yamato Rise and northwest 
in the coastal part of the Sea of Japan, the distance between 
the profiles is 10 miles. 

The survey was continued in 1991 and 1994 in the Central 
basin of the Sea of Japan (cruise No. 11 and 14 of the 
R/V “Professor Gagarinsky”). The observations were carried 
out by five GMN-K instruments with analog and then with 
digital recording. The error of observations in 1991 was 
+2.1 mGal, and in 1994 + 1.6 mGal. The direction of the 
profiles is meridional, the distance between the profiles on 
the 11th cruise was 10 miles, and on the 14th cruise —5 miles. 
Characteristically, since 1994, GPS receivers began to be 
used for navigational support for marine surveys of the POI 
FEB RAS, which significantly increased the accuracy of the 
survey. 

The Pervenets Rise is studied in detail in the 18th cruise 
of the R/V “Professor Gagarinsky” (1996). The survey was 
carried out by six GMN-K instruments with digital recording 
(error + 1.4 mGal). The direction of the profiles is northeast, 
the distance between the profiles is 5 miles. 

In 1996 (the 27th cruise of the R/V “Akademik M.A. 
Lavrentyev”), a single profile survey was taken of the north- 
ern end of the Central Basin and the Bogorov Ridge. The 
survey was carried out by six GMN-K instruments with 
digital recording (error + 1.4 mGal). 

The study of the northern part of the Central Basin 
continued on the 21st cruise of the R/V “Professor Gagarin- 
sky” (1997). The survey was carried out by five GMN-K 
instruments with digital registration (error + 2.0 mGal). The 
direction of the profiles is northwest, the distance between 
the profiles is 5 miles. 

In 1998, at the R/V “Professor Gagarinsky” (24th cruise), 
a regional gravimetric survey was carried out in the Central 
Basin using a series of meridional and latitudinal profiles. 
The distance between the profiles was 25 km. The observa- 
tion error was +1.6 mGal. 

In 2003, POI FEB RAS (cruise 36) of the R/V Professor 
Gagarinsky carried out a detailed geophysical (including 
gravimetric) survey of Peter the Great Bay. The average 
distance between the profiles did not exceed 2 km. The obser- 
vation error was + 1.4 mGal. The location of the observations 
was determined using GPS satellite system. 

In 2010, the R/V “Akademik M.A. Lavrentyev” (cruise 
52) of the POI FEB RAS carried out a gravimetric survey 
(42.0 mGal error) at the northern closing of the Central 
Basin, in the vicinity of the Vityaz Rise. 

This completed the first stage of areal geophysical work 
in the Central Basin of the Sea of Japan, within the eco- 
nomic zone of the Russian Federation. The root-mean-square 
error, calculated from the intersection of various gravimetric 
surveys, did not exceed 5 mGal. Measurement processing 
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was carried out according to standard methods. At the first 
stage, straight sections of the profiles were chosen, them the 
vessel moved with uniform speed and constant heading. The 
Eötvös correction, correction for the cross-dropping effect, 
for the inertia of the measuring systems of gravimeters, and 
zero drift were introduced. From the observed values, the 
normal gravitational field of the Earth, calculated according 
to the international formula of 1967, was subtracted. The 
anomalous field thus obtained was analyzed for intersections, 
after which gravimetric data catalogs were created. The 
construction of maps was carried out in the SURFER Golden 
Software program. The main operation at the initial stage of 
mapping was the creation of a regular grid of gravitational 
field values with a square cell (grid) based on an irregular 
pseudo-rectangular grid of observed data. When building the 
grid in the specified program, the so-called “search radius” 
was used, i.e. the distance within which data is searched 
to interpolate a regular grid value. In our calculations, the 
optimal search radius was set from 5 to 25 km depending 
on the degree of exploration of the water area, as a result of 
which three nearby profiles participated in the formation of 
the interpolation value of gravity at each grid node. The map 
constructed in this way is shown in Fig. 2a. 

Starting from 2017, the POI FEB RAS, at a new 
qualitative level, resumed geophysical exploration of the 
Sea of Japan. Marine gravimetric survey was carried out 
with Chekan-AM mobile gravimeters (Shelf-E modification) 
serial number 47. Studies were concentrated in the northern 
part of the Sea of Japan: the closure area of the deep-water 
Central Basin and the Tatar Strait. To date (October 2019), 
four expeditions have already been carried out in this area 
under the project “Integrated geological-geophysical, gas- 
geochemical and oceanographic studies in the Sea of Japan 
and the Tatar Strait.” 

So in 2017 at the R/V “Akademik Oparin” a geophysical 
survey was carried out in the southern part of the Tatar Strait 
from the latitude of the port of Sovetskaya Gavan in the north 
to the latitude of about. Moneron in the south. 112 geophysi- 
cal profiles with a total length of 3,860 miles were completed 
here. The survey error did not exceed +0.49 mGal. 

In 2018, two expeditions were completed. The research 
of the first expedition, to the R/V “Akademik M.A. Lavren- 
tyev”, was concentrated in the junction zone of the Central 
Basin of the Sea of Japan with the Tatar Strait trough valley. 
The second, at the R/V “Akademik Oparin”, passed north, 
from the latitude of the Soviet harbor to Cape Surkum. The 
survey error, in the first case, did not exceed £0.38 mGal, in 
the second +0.4 mGal. 

In May-June 2019, expeditionary studies were carried out 
in the remaining, not studied by us, part of the Tatar Strait, 
south of Moneron Island. During the expedition more than 
113 tacks were completed, with a total length of more than 
6,170 km. The survey error did not exceed +0.38 mGal. 
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Fig. 2 Map of gravitational 
anomalies in free air (conditional 
level), made according to ship 
on-board measurements for the 
period of expeditionary research 
from 1982 to 2010 (a) and from 
2017 to 2019 (b). The numbers 
indicate: 1 — Peter the Great 
Seamount; 2 — Bogorov Ridge; 
Rises: 3 — Pervenets, 4 — Vityaz, 
5 — Yamato, 6 — Alpatov, 7 — 
Lavrentiev 
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Gravimetric Studies in the Sea of Japan 
3 Results 


The works listed above made it possible to close white spots 
in the scientific geophysical study of the entire northwestern 
part of the Sea of Japan adjacent to the continent (Primorye, 
northern part of the Korea Peninsula). In recent years, the 
area of articulation of the Tatar Strait with the northern 
closure of the Sea of Japan and most of the Tatar Strait 
itself has been studied. This allowed, firstly, to build a series 
of conditional geophysical maps on a scale of 1:1,000,000 
(including gravimetric) and, secondly, to study the structural 
and genetic relationship of the Central Japanese basin with 
the margins of the adjacent continent. One of the directions 
of such studies was the gravitational (density) modeling of 
the earth’s crust at the junction of the indicated basin with 
the shelf and coastal geological structures of southern and 
south-western Primorye. 

Modeling is one of the effective methods for studying the 
deep structure, physical properties and processes occurring 
in the bowels of the tectonosphere, based on the use of 
available data and ideas about the environment. The authors, 
in particular, in the framework of these studies quite widely 
used the so-called gravitational (density, structural-density) 
modeling. A gravitational model is usually understood as a 
density model, the calculated gravitational field from which, 
up to a constant component, is adequate to the initial field 
of gravity (Bryansky 1991). “By calculating models of the 
density distribution in the geological environment that create 
an effect similar to the observed gravitational field, we can 
obtain information about the structure of the geological 
environment” (Krasovsky 1981). 

Structural-density modeling (Fig. 3) was carried out 
according to a series of profiles crossing the transition 
zone from the edge of the continent to the central deep- 
sea basin of the Sea of Japan, which quite clearly showed 
the characteristic features of the transformation of the 
continental crust into the oceanic crust. 

Transformation of the crust begins on the continent at a 
considerable distance from the upper edge of the continental 
slope. In different sections, this distance is not the same. 

The transformation of the continental crust is both struc- 
tural and substantial nature. The main factor determining 
this transformation is the reduction in the thickness of the 
“basalt” layer from the continent towards the basin of the 
Sea of Japan and its replacement with a mantle substrate. 
The supra-basaltic crust (crystalline basement) in most cases 
practically does not change its thickness to the continental 
slope. However, the substantial composition of this founda- 
tion in different parts of the studied area does not remain the 
same. 

The final transformation of the continental crust occurs 
within a narrow strip of the outer shelf — the bottom of 
the continental slope. Here its sialic part disappears com- 
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pletely, the thickness of the lower crust (“basalt” substrate) 
is significantly reduced, which is replaced by mantle masses 
from below. In all the sections presented in the slope part, 
an increase in the thickness of the “transition” layer is 
observed, which indicates the processes of effusion of effu- 
sive rocks that took place here, with the formation, in some 
cases, clearly expressed in the relief of volcanic structures 
(Seamount Peter the Great). 

At the base of the slope and the slope of the deep-sea 
basin, a narrow zone of deformations of all layers of the 
oceanic crust is recorded on all profiles, which can be a 
morphological expression of a powerful tectonic structure 
(tectonic seam) at the junction of dissimilar types of the 
earth’s crust. 

Another area of research was the depth’s determination 
of the Mohorovičić discontinuity (Moho surface) bedding 
within the Sea of Japan. It is well known that on the 
Mohorovičić discontinuity there is a sharp increase in the 
speed of seismic waves and, accordingly, the density of the 
underlying rocks. Similar changes in density are observed 
during the transition from the water column to the upper part 
of the sedimentary layer and during the transition from the 
lower part of the sedimentary layer to consolidated rocks of 
the underlying layers (acoustic foundation). The observed 
gravitational field reflects the integral effect of the above 
boundaries and can be used to determine the depth of their 
occurrence by statistical methods. To determine the depth of 
the Mohorovičić discontinuity within the study area, we used 
the statistical dependence of the magnitude of the averaged 
gravitational anomalies, sedimentary cover thickness and 
sea depth on the depth of the Mohorovičić discontinuity, 
proposed in (Su Datsuan 1982): 


H = 33.49 — 0.063Ag free — 0.00482 Hsea — 0.001 7A sea 


where H is the depth to the Mohorovičić discontinuity, km; 
A free — anomaly in free air, mGal; Hsea — depth of the 
seabed, m; Asea — thickness of the sedimentary layer, m. The 
depths of the seabed (Hsea) were determined according to the 
sonar data, seq were determined according on the NSP data. 

As is known, there is a statistical relationship between 
the thickness (the depth of Mohorovičić discontinuity) and 
the type of the earth’s crust (Belousov and Pavlenkova 
1985). Moho depths in the indicated area are shown in the 
Figure 4. As can be seen from the figure, the relief of the 
base of the crust of the Central Basin of the Sea of Japan 
and its neast framing is of considerable complexity, which 
indicates the heterogeneity of different parts of this region. 
According to this feature, the studied area can be divided, 
first of all, into two large sections, the border between which 
runs approximately along the meridian of the underwater 
Pervenets Rise (132°30/). To the east of this boundary is 
the deepest part of the basin with a maximum protrusion of 
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Fig. 3 Structural-density models of the earth’ crust in the transition formations, 8 — granites, 9 — faults, 10 — a) observed gravitational field, 
zone from the mainland to the Sea of Japan. 1 — water column; 2 — b) calculated gravitational field, c) magnetic field. Numbers — densities 
sedimentary layer; 3 — volcanic sedimentary layer; 4 — “basalt” layer; (g/cm?) of rocks. The inset shows the position of the profiles. RF — 
5 — crystalline basement, 6 — gabbroids, 7 — pre-Mesozoic sedimentary Russian Federation, PRC — People’s Republic of China 
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Fig. 4 The map of the Moho discontinuity isodepths (kilometers below 
the sea level) and the types of crust beneath the north-western part Sea 
of Japan based on the marine gravimetric data. 1 — continental crust; 


the Mohorovičić discontinuity, the depth of which from the 
day surface varies from 14 km in the west to 12 km or less 
in the east. Without a water layer, the crust thickness here 
is 10.5-8.5 km, respectively, which is in good agreement 
with the results of seismic studies (Hirata et al. 1992). The 
Earth crust here is of the oceanic type. This area has a 
wedge-shaped east-north-east strike, complicated by local 
thickening of the crust up to 16—20 km within the underwater 
elevations. 

The western half of the basin differs significantly from the 
previous one by the great depths of the Mohorovičić discon- 
tinuity (16-18 km) and a different pattern of its relief. The 
central place here is occupied by a vast, laid out section of the 
western end of the basin of almost isometric shape, within 
which the depth of Moho is 15—16 km. Two “apophyses ” 
with a relatively high position of Moho, separated by an East 
Korean Rise, depart from this section in the south-west and 
south directions. The thickness of the consolidated crust here 
is 13—15 km, which allows us to attribute it to the suboceanic 
type. 

The thickness of the crust within large submarine hills 
reaches 26 km, which is characteristic of the Earth’s crust 
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2 — reduced continental crust; 3 — oceanic crust; 4 — presumable oceanic 
crust with the thickened sedimentary cover (suboceanic?) 


of a subcontinental type. The transition from the deep- 
sea basin to the continent is accompanied by the intense 
lowering of Moho to depths of 30 km and more, which 
corresponds to the minimum thickness of the continental- 
type crust. 


4 Conclusion 


For more than half a century, Russian and foreign researchers 
have been studying the geological structure and natural fields 
of the Earth in the waters of the Sea of Japan. Scientists 
of the Pacific Oceanological Institute of the Far Eastern 
Branch of the Russian Academy of Sciences made the most 
complete contribution to this process, in the waters of the 
economic zone of Russia (formerly the USSR). Gravimetric 
studies, as one of the deepest geophysical methods, made 
it possible to determine the thickness and types of the 
earth’s crust, study its internal substantial-density structure, 
separate blocks and trace tectonic faults. The integration of 
gravimetry with other geophysical methods, especially with 
deep seismic sounding, significantly improved the accuracy 
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of these studies and reduce the ambiguity of the result. Thus, 
we can say that the first stage of the study of the gravitational 
field of the Central Basin of the Sea of Japan has been 
completed. In the future, it is necessary to focus on a detailed 
study of individual morphostructures, the implementation of 
near sea bottom gravimetric observations and the gravity 
field’s gradients calculation, which, in combination with 
other geophysical methods, will provide new information on 
the structure and substantial composition of the earth’s crust 
in the region. 

The fundamental study and interpretation of the gravita- 
tional field in the waters of the Tatar Strait, in the northern 
part of the Sea of Japan, has just begun. The four expeditions 
carried out, unfortunately, did not allow us to characterize 
the gravitational field of this morphological structure with 
the required detail, especially in the area of transition from 
the Central Basin of the Sea of Japan to the southern part of 
the Tatar Strait. The authors hope that research in this area 
will be continued, including in cooperation with Japanese 
colleagues. 
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Abstract 


The article presents the description of two algorithms used for processing of the raw data 
of a gravity gradiometer. These algorithms are intended for estimation of some instrument 
errors. The first algorithm is applicable for the instrument operation in its stationary mode, 
the second proposes the use of a special test bench. Rotary gravity gradiometer of the 
accelerometric type was taken as a prototype for relevant mathematical models. Nowadays 
this type of gradiometer is brought to the stage of practical implementation and serial 


production. 


Keywords 
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1 Introduction 


The gravity gradiometer, that is an analogue of the model 
developed by Bell Aerospace Laboratory (USA) and cur- 
rently produced by Lockheed Martin (USA), is considered 
(Dransfield et al. 2010; Murphy 2010; Hammond and Mur- 
phy 2003). 

The instrument design is based on the slowly rotating 
disk with diameter equals 0.2 m and angular rate 0.25 
Hz. The parallel pairs of accelerometers are set on the 
diametrically opposite sides of the mentioned rotating disk 
and at the same distance from the center of the disk. The 
accelerometer sensitivity axes are tangent-oriented (Fig. 1). 
Further, for certainty, we assume that disk has the vertical 
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axis of rotation. This construction of a gradiometer is called 
GGI (gravity gradiometer instrument). 

Gradiometer designed for airborne gravimetry. Gradiome- 
ter should be installed on a gyro stabilized platform. 

The linear combination of accelerometer readings gener- 
ates the GGI output signal. The rotation of the disk induces 
forced harmonic oscillations, and as a result the output signal 
becomes modulated at a frequency equal to twice angular 
rate: 


W=(fith)-(B+ fy = 
= I[(Px2 — Tii + o2 — œ) sin 20t+ (1) 
+ (2 + ww) cos 221). 


Here 
— W is the GGI output; 
— fi,i = 1,...,4 are the components of a specific force 


measured by accelerometers 1, 2, 3, 4; 

— Di, I», Tiz are the components of the gravity gradient 
tensor l’ = grad g, where g is the vector of the specific 
gravity; 

— 2 is the angular rate vector of the disk, and 22 > u, 
where u is the absolute velocity of earth’s rotation; 

— a is the absolute angular rate vector of a gyro platform 
(GGI body frame). 
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Spin axis 


Asselerometer 


Rotating disk 


Fig. 1 GGI model 


The problem of estimation of J” tensor component using 
W measurements solved in post-processing. In this case one 
have to use current navigation parameters of high accuracy, 
in particular w. 

The basic method used for increasing the accuracy of the 
GGI operation is the compensation of its instrument errors. 
The most stable parameters of these errors can be estimated 
during calibration that proposes the usage of a specialized 
test bench. There were discussed several options of the 
mentioned calibration procedure in recent publications. The 
first variant proposes the usage of a centrifuge (Yu abd Cai 
2018). The second variant is based on the special motion of 
the disturbing mass (Deng et al. 2018). 

The pre-start calibration mode of GGI is also possible, 
that consists of the regular operation at a fixed point for a 
short time. 

It should be noted that the gradiometer error model 
includes the combinations of the instrumental errors of the 
paired accelerometers. So, it is important to do estimation of 
these combinations exclusively, and not to do estimation of 
errors of single accelerometers. This is why the calibration 
procedure should be performed for the assembled GGI. 


2 Model of Instrumental Errors of the 
Gradiometer 


The accepted model of instrumental errors is as follows. Let 
f be the measured value, f’ be the result of the measure- 
ment, then f’ = f + Af, where Af is the instrumental 
error. We will use the following designations (Golovan et al. 
2018): 
— & = (61, 62,83)" is the small rotation vector that charac- 
terizes the installation errors of a GGI on a gyro platform; 
— Al; = (l; —1)/2 is the inaccuracy in the distance from 
the center of the proof mass of the i -th accelerometer with 
respect to the center of the rotating disk (the true distance 
equals to /; /2), 
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— Xi is the misalignment error of sensitivity axis of the i-th 
accelerometer; 
— A fio is the bias of the i-th accelerometer measurement; 
— ki; is the scale factor error of i-th accelerometer measure- 
ment; 
— Afis is the noise component of accelerometer measure- 
ment. 
The error model of the i-th accelerometer is as follows 
ff =U+k) fi t+ Afio+ Afis + xi fis (2) 
where the component y; fi reflects the cross-coupling 
effect of the specific force acting on the proof mass of 
the accelerometer orthogonally to its longitudinal axis due to 
a skew in the disk’s plane. 
Let W” = (f/ + f2) — (J3 + fi) be the measurement of 
the W (1). Measurement error AW = W’ — W contains the 
frequency components Ag, A1, A: 


AW = Aj + A; + Ad. (3) 

Taking into account the mentioned definitions, the error 
model of the gradiometer output takes the form (Golovan 
et al. 2018): 


Ao = (A fio + A foo — A foo — A fao)+ 


+ (A fis + A fos — A fos — A fas)— (4) 


1, 
— 5@3| ki + ka — ka — ka + xı + x2- x3 — at 


4 Ah + Al, — Als — =, 
l 
Ay = | (ki — ko + %3 — X4)w2+ 
+ (kz — k4 — x1 + x2)wı | cos 2t— 
= [k — k + %3 — X4)wı— 
— (k3 — k4 + Xi — x2)w2] sin Èt, 


(5) 


where w1,w2 are the horizontal components of a gyro plat- 
form’s specific force, © is the vertical component of the 
angular acceleration, 
Ay = 2l [61 (I3 + w3) + ô&2(Ti3 + @103)— 
— 263(1\2 + @1>)| sin2Qt+ 
+ 21 [283 (M2 — Ty +0? — 07) + 83 + @1@3)— 
— 69(193 + 23) | cos 22t+ 


+ (ki the -ks ka + %1 + X2- X3 — Xa+ (6) 


Al Ah — Al — Al 
42 it ~ 3 *) x 
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l 
x RE -Õli + ws — w?) sin 22t+ 


l 
+ z702 + ww) cos 2021]. 


3 Calibration Algorithm 


The following combinations of the GGI instrument error are 
of interest: 

yı = ĝi, y2 = ĝ2, y3 = 63, 

y4 = A fio + A fro — A fo — A fao, 

p=k tkk ktat- | 


Ahi + Al, — Ahl — Aly 


It is obvious that the constructive instrument errors — 7; 
(that are the mutual misalignments of the sensitive axes of 
accelerometers) and Al; (that is the inaccuracy of disk’s 
radius value) remain unchanged in GGI operation, so it is 
necessary to include them in the relevant calibration problem. 
At the same time ô; (that is the GGI installation error), A fio 
and k; (that are the biases and scale factor errors) may vary 
from run to run. So these parameters should be estimated 
during GGI operation. 


3.1 Case of Fixed Point 

The following reference information is used for calibration 

of the gradiometer at a fixed point: 

— the vector of earth rotation with respect to the geographi- 
cal reference frame 


OQ, = 0 
w2 = UCOS@ (8) 
wz = using 
— the vertical component of angular acceleration equals 
Zero: @3 = 0; 
— the linear accelerations equal zero w; = 0, w2 = 0; 
— the reference values for the gravity tensor components 
Ty. 
Then the model for output error of the gradiometer takes the 
form 


z=W-W=AW=Al(A+ fA)-(A+ fl = 
= y4 F (A fis Ea A fas a A fs = A fas)— 


+ 21| y1 (r + u? sing coso) + y2I\3 — 2y3Ti2+ 
l 
+ gs02 — I, +u’ cos? 0] sin2Qt+ (9) 
+ 21 [273z — Dı +u? cos? o) + yi Ti3— 
2. l 
— yo(Io3 + u sing coso) + syst cos29t. 


Denote by x1, x2, x3 the amplitudes of the harmonics in 
the z measurement: 
x)= y4 (10) 


xXx = 2l [nis + u? sing cosg) + y2I3— 2vsM2|+ 


I 
+7rs 2 — Du + u’ cos” g) (11) 
x3 = 2l [2732 — Tii + u? cos” o) + yıli3— 
20: l 
—y2(I3 + u sing cosg) | + 5 Yt ia- (12) 


Let assume that instrumental errors are constant during the 
experiment: 


Afi =0,i =1,...,4 
§ =0,i=1,...,3 


k; =0,i=1,...,4 (13) 
4% =0,i=1,...,4 
Ai; =0,i=1,...,4. 


Then the estimation problem for the state vector x = 
(x1, X2, x3)” of linear dynamic system can be set as 


x = 0, 


II 


z = xı +x. sin 22t + x3cos22t +r, (14) 
where r is the white noise of a given intensity that is 
composed by the linear combination of noises of single 


accelerometers 


r= A fis Æ A fs a A fis == A fas (15) 
The problem posed can be solved by application of the 
relevant Kalman filter. 
Since the set of functions 1, sin2 2t, cos 2S2t are linearly 
independent, the constant factors in these functions (14) 
become observable. 
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So, the estimates of combinations (7) can be obtained 
and above algorithm can be used in as a pre-start calibration 
algorithm. 


3.2 The Usage of a Stand of Linear 


Movement 


The estimates for combinations kı + k2 —k3—k4 + 741+ X2— 
X3 — X4 can be derived using test bench of linear movement. 
Calibration procedure proposes that the gradiometer moves 
in the horizontal plane with known linear accelerations along 
a given direction. The algorithm does not require reference 
values of the gravitational gradient tensor components in this 
case. 
One can form measurements in the following way 


Z = fi- h = A fio — A fro E A fis = A fos— 


= [(k1 + k2)w, + Gal + X2)W2] sin 2Qt+ (16) 
+ [(kı + k2)w2 + (x1 + X2)wi] cos 2t 

z = fh — fi = Afso —A fy + A fas — A fas— 
— [3 + k4)w2 — (13 + x4)wı] sin 2t— (17) 


— [(k3 + k4)wı + (73 + X4)w2] cos 2t. 


When one selects the orientation of the stands rail in 
the North-South direction (the linear acceleration component 
w, = 0), then the measurements zı, z2 take the form: 


H = A fio = A fro + A fis _ A fos— 
— (1 + X2)w2 sin Zt + (kı + k2)w2 cos Qt, 


(18) 
2= A fro T A fao ae A fas = A fas— 
— (k3 + ka)w2 sin Qt — (%3 + X4)w2 cos Lt. 
Let introduce: 
xı =A fio — A fro x3 = A fzo — A fao 
xo =k, + ko X4=k3+ k4 (19) 
x3 = X1 + x2 X5 = X3 + X4. 
Then 
zı = X1 + x2w2 cos Mt — x3w2 sin Üt + rı 
(20) 


Z = X4 — X5w2 Sin Rt — xXew2 cos Lt + r2. 


Thus, the problem is reduced to two independent estima- 
tion problems 


x; =0 x4 = 0 
x2 = 0 and x5 =0 (21) 
x3 = 0 xX =0 


with the help of aiding measurement z1, z2. 
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Since the functions 1, sin 2t,cos Mt are linearly inde- 
pendent, then the constant factors in these functions are 
observed. Therefore, the following combination kı +k2, k3 + 
k4, X1 + X2, X3 + X4 can be estimated, e.g. applying Kalman 
Filter. 

The similar algorithm can be implemented for the estima- 
tion of kı + kz and k3 + k4 combinations in GGI regular 
operation mode (O’ Keefe et al. 1997). 

The y1, y2, Y3, Y4, ¥5 Components (7) in the GGI error 
model (4, 6) can be compensated on the basis of correspond- 
ing estimates derived by Kalman filters. 


4 Conclusion 


The algorithms were elaborated for the estimation of instru- 
ment errors of GGI. The implementation of these algorithms 
is fairly easy, they can be included into calibration procedure. 
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Abstract 


In development of a rotating accelerometer gravity gradiometer (RAGG), it is difficult for us 
to distinguish the measurement signals and error components in the RAGG output without a 
verified and correctness RAGG analytical model. In addition, many key techniques, such as 
RAGG analytical model validation, online error compensation, post motion error compen- 
sation, are difficult to be verified. RAGG numerical model can provide validation platform 
for solving all the above problems, which can speed up the development of the RAGG. 
In this study, based on the principle and configuration of the RAGG, we synthetically 
consider almost all the error factors, such as circuit gain mismatch, installation errors, 
accelerometer scale-factor imbalance, and accelerometer second-order error coefficients, 
construct a parameters adjustable RAGG numerical model. In multi-frequency gravitational 
gradient simulation experiment, we use the RAGG numerical model simulating the situation 
that a test mass rotates about the RAGG with time-varying angular velocity to generate 
multi-frequency gravitational gradient excitations; the experiment results are consistent with 
the theoretical ones; the RAGG numerical model can recur some phenomenons of a actual 
RAGG. 


Rotating accelerometer gravity gradiometer - Numerical model - Center gravitational 
gradients - Virtual rotating accelerometer gravity gradiometer 


superconducting gravity gradiometers, cold atomic interfer- 
1 Introduction ometer gravity gradiometers (Paik 2007; Difrancesco 2007; 
Moody 2011; Liu 2014; Hao 2013). Rotating accelerometer 
Airborne gravity gradiometry is an advanced technology gravity gradiometer (RAGG) was developed by Ernest Met- 
for surveying a gravity field; it acquires gravity field in- zger of Bell-Aerospace in the 1980s (Heard 1988). But, over 
formation with high efficiency and high spatial resolution. the past decade, rotating accelerometer gravity gradiometers 
Compared with gravity information, the gravity gradient are the only commercial moving base gravity gradiometer 
tensor provides more information on the field source such in the word; all other types are either in fight testing or 
as orientation, depth, and shape (Tang and Hu 2018; Yan in a laboratory setting. Companies that operate commer- 
and Ma 2015). There are many different types of gravity cial rotating accelerometer gravity gradiometer systems are: 
gradiometer: rotating accelerometer gravity gradiometers, Bell Geospace (Air-FTG), ARKeX (FTGeX), GEDEX (HD- 
AGG), and FUGRO (AGG-Falcon) (Rogers 2009; Metzger 
1977). In this study, we synthetically consider almost all 


M. Yu - T. Cai (è<) a a the unperfect factors, for example accelerometer installa- 
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Nanji i tion errors, accelerometer scale-factor imbalance, etc., and 
anjing, China : ; 

e-mail: caitij@seu.edu.cn construct a numerical model. The RAGG numerical model 
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is a virtual RAGG with a comprehensive set of precisely 
adjustable parameters; based on it, many key techniques, for 
example RAGG automatic calibration, post error compensa- 
tion, and self-gradient modeling, can be verified. 


2 RAGG Numerical Model 


A high-precision numerical model of the RAGG is estab- 
lished, as shown in Fig. 1. In the numerical model, each 
accelerometer has six mounting error parameters: radial dis- 
tance, initial phase angle, altitude angle, and misalignment 
error angles. Among them, the radial distance, initial phase 
angle, and altitude angle determine the mounting position 
of the accelerometer; the misalignment error angles deter- 
mine the orientation deviation between the accelerometer 
measurement frame and the accelerometer nominal frame 
of the actual mounting position. The detailed definition of 
the mounting error parameters can be found in literature (Yu 
and Cai 2019). Moreover, each accelerometer has nine other 
output model parameters: zero bias (K jo), linear scale factors 
(Kj 1), second-order error coefficients (K j2, Kj4, Kjs, Kjo, 
Kj7, Kj), and current to voltage gain (k jy 7). We use a test 
mass to produce gravitational gradients to excite the RAGG. 
The specific force of the accelerometer A; in the RAGG 
numerical model is given by: 


fi = Semm + @im X Pon Aj + @im X (@im x vom A;) 
ADD 
-GmA;S f|A;S} . 


Where f emm is the specific force of the RAGG; @;,, is the 
angular acceleration of the RAGG with respect to the inertial 
frame; @;m is the angular velocity of the RAGG with respect 
to the inertial frame; G is the gravitational constant; ro, 4; 
is the position vector of accelerometer A; in the RAGG 
measurement frame; m represents the weight of the test 
mass; and A; S is the position vector from accelerometer 
A; to the test mass. If the test mass is not a point mass, 
the gravitational acceleration that the RAGG accelerometers 
undergo produced by the test mass can be calculated using 
finite element analysis. In addition, the test mass can be 
in motion with respect to the RAGG; in this case, A; S is 
time varying. The specific forces of accelerometer A; in the 
accelerometer nominal frame of the actual mounting position 
(fix. fiy» jz) can be calculated from: 


fix = fj Tix» 
fiy = Sj Tiy 
Fic = SF Tic. 


(2) 
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sin (2Qt) 


REF Signal 
COS (2.21) 


QAM 
Demodulator 


Inline channel +> 


BPF Cross channel > 


1. setting RAGG simulation parameters 
a. accelerometer mounting parameters 
b. accelerometer model parameters 

c. RAGG motion parameters 


(a) 


d. simulation duration time 


calculate the specific force f; f, f, in accelerometer 
measurement frame at time t 


calculate accelerometer voltage output at time t 


calculate RAGG output before demodulation 
at time £ Gou(Q=Vi(O+V(0-V3(0-VaD 


Is time <= duration? 


QAM demodulator 


RAGG output after demodulation 


(b) 


Fig. 1 Principle and program flow of the RAGG numerical model. 
(a) Principle of the RAGG numerical model. (b) Program flow of the 
RAGG numerical model 
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Where T jx, T jy, and T j; are unit vectors of the accelerom- 
eter nominal frame of the actual mounting position in the 
directions of the x-, y-, and z-axes. The specific forces in the 
accelerometer measurement frame are: 


fii fix 
fio| =C | fiy |. (3) 
fiv fiz 


where C is the transformation matrix from the accelerom- 
eter nominal frame of the actual mounting position to the 
accelerometer measurement frame. To make the numerical 
model approximate the actual RAGG, we add accelerometer 
noise to the accelerometer model: 

Vie 2 2 
Emaka Íinoise + fji + Kjo + Kyo fii + Kjsfjo 
+K; fip + Kishi fjo + Kishi fip + Kis fiofip - 

(4) 


The accelerometer noise fjnoise is simulated by a power 
spectral density model (Jekeli 2006): 


D(f noise = af + Or, (5) 


where œ and b represent the amplitude and low-frequency 
growth of the red noise, and wr denotes the amplitude of the 
white noise. 

Figure |b is the program flow of the RAGG numerical 
model. Firstly, the RAGG simulation parameters are set up, 
including test masses parameters, RAGG rotating disk pa- 
rameters, accelerometer mounting parameters, accelerometer 
model parameters, RAGG motion parameters, etc. Then sub- 
stituting the parameters into the formula (1)—(3) calculates 
the specific force in accelerometer measurement frame at 
time ¢. According to the formula (4), calculating the out- 
put voltage of the RAGG accelerometer, the RAGG output 
before demodulation at time t is calculated by: Gou (t) = 
V(t) + Vo(t) — V3(t) — Va(t). The above process is repeated 
until time ¢ is equal to the simulation duration time. Finally, 
the RAGG output data is input to the quadrature amplitude 
modulation (QAM) demodulator to extract gravitational gra- 
dient. 

Let Ix, Ixy, Iz, My, and Ty: represent the five inde- 
pendent gravitational gradient elements at the origin of the 
RAGG measurement frame. When mass is far enough away 
from the RAGG, the gravitational acceleration measured 
by the RAGG accelerometers is a first-order approximation 
of the gravitational acceleration and gravitational gradient 
tensor at the center of the rotating disc; in this case, the 
inline channel measurement and the cross channel mea- 
surement of the RAGG approximate I, — I, and Iyy; 
otherwise, the inline channel measurement and cross channel 
measurement of the RAGG are the sum of I, — Iyy, 
Ty, and high-order gravitational gradient tensor elements. 


To distinguish between I, — Ixy, [xy and the measure- 
ments of the RAGG, I, — Ixy, I, is called as center 
gravitational gradients; lyx — I, is the inline channel of 
the center gravitational gradients; ly is the cross channel 
of the center gravitational gradients. In RAGG analytical 
model, gravitational acceleration that the RAGG accelerom- 
eter undergoes is a first-order taylor approximation of the 
gravitational acceleration and gravitational gradient tensor 
at the center of the rotating disc; but in the numerical 
model, the gravitational accelerations are calculated using 
Newton’s law of gravitation instead of a linear approxima- 
tion. Therefore the numerical model are close to the actual 
RAGG. 


3 Multi-frequency Gravitational 
Gradient Simulation Experiment 


In multi-frequency gravitational gradient simulation 
experiment, a test mass rotates about the RAGG with 
time-varying angular velocity producing multi-frequency 
gravitational gradient excitations. Based on the angular 
velocity of the test mass and its initial coordinate in the 
RAGG measurement frame, we can obtain the coordinates 
of the test mass in the RAGG measurement frame at 
any time. We then calculate the gravitational gradient 
tensor at the origin of the RAGG measurement frame 
and then calculate the center gravitational gradients. 
The RAGG numerical model simulates a perfect RAGG 
without accelerometer mounting errors, accelerometer 
scale-factor imbalances, accelerometer second-order 
error coefficients and accelerometer noise, so we set 
the accelerometer mounting errors, the accelerometer 
second-order error coefficients, accelerometer noise 
parameters to zero. The linear 
four accelerometers is kit = 10 mA/g, the current- 
to-voltage gain is kjyjp = 10° ohm, the nominal 
mounting radius R is 0.1 m, and the rotation frequency 


scale factor of the 


Fig. 2 A test mass rotating about the RAGG with time-varying angular 
velocity 
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Fig. 3 Accelerometer voltage outputs in RAGG numerical model 
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Fig. 4 Output voltage of the numerical model before demodulation 


of the RAGG disc is 0.25 Hz. Figure 2 shows a point 
mass of 486 kg with an initial position in the RAGG 
measurement frame of (1.1,0.5,0) and rotating about 
the RAGG with time-varying angular speed w(t) = 
500 + 400 sin(0.0628r)° /h. 

Figure 3 shows the voltage outputs of the four accelerom- 
eters in the RAGG numerical model excited by the rotating 
point mass. Figure 4 shows the output voltage before demod- 
ulation of the numerical model. Figures 5 and 6 show the 
demodulated gravitational gradient comparison among the 
RAGG model and the center gravitational gradients; we can 
see that the inline channel and the cross channel of the RAGG 
numerical model are overlapped with that of the center 
gravitational gradients and only one curve is displayed; to 
distinguish the outputs of numerical model and the center 
gravitational gradients, their differences are calculated and 
shown in Figs.7 and 8; the differences are in the order of 
10 7! E, and they are caused by the high-order gravitational 
gradient tensor. 
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Fig. 5 Demodulated gravitational gradient comparison between nu- 
merical model and center gravitational gradients: Inline channel 
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Fig. 6 Demodulated gravitational gradient comparison between nu- 
merical model and center gravitational gradients: Cross-channel 
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Fig. 7 The difference between center gravitational gradients and nu- 
merical model: Inline channel 
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Fig. 8 The difference between center gravitational gradients and nu- 
merical model: Cross-channel 


4 Conclusion 


Based on the measurement principle and configuration of the 
RAGG, we considered the factors of circuit gain mismatch, 
installation error, accelerometer scale-factor imbalance, and 
accelerometer second-order error coefficients, then devel- 
oped a high-precision numerical model. The parameters of 
the RAGG prototypes are unknowable and uncontrollable, 
but the parameters of the RAGG numerical model are ad- 
justable and knowable; to some extent, the RAGG numerical 
model are more suitable for verifying some technique solu- 
tions of the RAGG. 
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Abstract 


The paper presents an approach to determination of the gravity disturbance vector from 
airborne gravimetry measurements at an aircraft’s flight path. A navigation-grade inertial 
navigation system (INS) and the carrier-phase differential mode of GNSS are assumed. 
To improve observability of the gravity horizontal components, which are observed in 
combination with the INS systematic errors, we use a spatial model of gravity. We 
parameterize the disturbing potential in the observation area using the spherical scaling 
functions. The unknown coefficients of the parameterization and the INS systematic 
errors are estimated simultaneously using the Kalman filter. Due to ill-conditioning of 
the estimation problem, the information form of the Kalman filter and regularization 
are used. The numerical results obtained from simulated data processing shows that the 
approach based on spatial modeling is capable to improve accuracy of the gravity horizontal 


component determination comparing to a typical modeling of gravity in the time domain. 


Keywords 


Airborne gravimetry - Gravity disturbance vector - 


functions 


1 Introduction 


Airborne gravimetry is capable to provide Earth’s grav- 
ity observations of high accuracy and spatial resolution. 
A typical airborne gravimetry system includes an inertial 
navigation system (INS) (strapdown or a gyro-stabilized 
platform) and global navigation satellite system (GNSS) 
receivers operating in the carrier-phase differential mode. 
Airborne gravimetry is well established for determining the 
vertical component of the gravity disturbance vector (scalar 
gravimetry based on platformed INS) and is capable to 
provide all three gravity vector components simultaneously 
from airborne measurements (vector gravimetry) (Kwon and 
Jekeli 2001; Becker et al. 2016). The common method used 
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Kalman filter - Spherical scaling 


in vector gravimetry is based on INS/GNSS integration 
via Kalman filtering. A high accuracy (navigation grade) 
strapdown INS is assumed. The main difficulty is in deter- 
mining the gravity horizontal components (deflections of 
the vertical) as these are observed in combination with 
the INS systematic errors (the inertial sensor biases, the 
attitude errors, etc.), which have a similar spectral content. 
An additional information (e.g., an a priori gravity model, 
external gravity data) is required for separation. 

The common approach to airborne vector gravimetry is 
based on modeling gravity in the time domain during Kalman 
filtering, e.g., assuming each gravity component to be a 
stationary stochastic process (Becker et al. 2016; Jensen et al. 
2019). As an alternative, a gravity model may not be used 
explicitly, and the gravity components can be determined 
from the Kalman filter residuals (Kwon and Jekeli 2001). 
However, remains of the systematic errors may be present 
in the gravity estimates of both approaches, and additional 
corrections should be applied. 
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In this study, we develop an approach that takes into 
account spatial correlation of gravity at adjacent survey lines 
to improve observability of the gravity horizontal compo- 
nents. The approach is based on spatial modeling of the dis- 
turbing potential in the observation area using the spherical 
scaling functions (SSF) of a certain resolution level, which 
are harmonic and spatially localized (decrease with distance 
from their origins) (Freeden and Michel 2004). The SSFs 
and other spherical radial basis functions are widely used 
in local and regional gravity field modeling and are mainly 
applied for processing satellite gravity data and combining 
gravity data, e.g., Lieb et al. (2016) and Klees et al. (2008). 
We use the SSFs to airborne vector gravimetry problem 
solving. The developed estimation algorithm is based on 
simultaneous estimation of the INS systematic errors and the 
scaling coefficients (the coefficients of the parameterization 
of the disturbing potential) via Kalman filtering given mea- 
surements at the aircraft’s flight path. The information form 
of the Kalman filter (Kailath et al. 2000) and regularization 
are applied due to the ill-conditioning of the estimation 
problem. The approach was partially presented in Bolotin 
and Vyazmin (2016). 

The aim of this work is to investigate the performance 
of the developed approach by processing simulated airborne 
data and by comparing with the common approach based on 
modeling the gravity vector components in the time domain. 
The numerical results obtained with the developed approach 
are optimistic and show that an increase in the accuracy of the 
gravity horizontal component determination can be achieved. 


2 Airborne Vector Gravimetry 
Equations 


The gravity disturbance vector determination is based on 
INS/GNSS integration. The basic equation is the equation of 
motion of the INS proof mass M expressed in the navigation 
(local-level) frame M x as 

dy = OF, +Â + LE fet ystig D 
Here v, is the 3 x 1 vector of the aircraft velocity relative 
to the Earth expressed in Mx frame coordinates, ôg, 
is the gravity disturbance vector, y, is the reference 
(normal) gravity vector, f, is the specific force expressed 
in coordinates of the body frame (denoted by M z). Further, 
Lx is the transformation matrix from the navigation frame 
(Mx) to the body frame which satisfies the kinematic 
equation (Farrell 2008), Ors is the (absolute) angular 
velocity of the navigation frame with respect to the 
Earth-centered inertial frame, expressed in Mx, OnE 
is the Earth absolute angular velocity vector expressed 
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in Mx, One is a skew-symmetric matrix formed by 
the components of the vector Ore in such a way that 
Otx = 
product. 

From Eq. (1), we arrive at the INS mechanization equa- 
tions (Farrell 2008) by replacing true values in Eq. (1) with 
measurements from the INS sensors (measurements of the 
body-frame angular velocity and the specific force) and 
with GNSS position and velocity, and omitting the gravity 
disturbance vector term. Here, we assume GNSS positioning 
and velocity solutions to be obtained using the carrier- 
phase differential mode. In this case, inaccuracies in GNSS 
positions can be neglected (Vyazmin and Bolotin 2019). 
Denote by v’, the velocity solution obtained from solving the 
INS mechanization equations and by L/, the transformation 
matrix from the navigation frame to the computed body 
frame. 

By subtracting Eq. (1) from the INS mechanization equa- 
tions, the INS error dynamics equations can be obtained as 
follows (Farrell 2008; Vyazmin and Bolotin 2019) 


Ors xX vy, Where x means the cross 


+ Li Af - ò’ L'AO, 
Ax = Oa, — Li Aw (3) 


where œ, is the attitude error, dv, is the velocity error 
defined as the sum of the total velocity error Av, = 
v — vx and the term 0,a,. Note that the equations do 
not contain the INS positioning solution as we assume that 
the INS proof mass position is perfectly determined from 
GNSS. 

The terms Aw, Af are the gyroscope and accelerometer 
errors, respectively, for which the following models are 
assumed in the study: 


(4) 
(5) 


Aw = ba + 4%, Af =bf+qpr, 


ARNE NETA 


where bo, by are the gyroscope and accelerometer biases, 
respectively, the vectors 4o, 4f, €, €f are the random 
errors. The error models in Eq. (5) represent the sensor bias 
drifts. 

Denote by v5 the GNSS velocity of the aircraft 
relative to the Earth and by y, the observation computed as 
y, = v’, — vo%SS_ The observation equation is written as 
(Vyazmin and Bolotin 2019) 

(6) 


al 
Yx = Vx — Vay — ey 


where e, is the GNSS velocity error. 
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3 Local Gravity Model and Estimation 
Algorithm 


In order to use a Kalman filter framework, we introduce 
a model of the gravity disturbance vector. We use a local 
spatial model by parameterizing the disturbing potential at 
the points at the aircraft’s flight path by the SSFs (Freeden 
and Michel 2004; Klees et al. 2008) 


N 
T(r)= X aidr) 


i=l 


(7) 


where a;, i = 1,..., N, are the unknown parameters (the 
scaling coefficients) to be estimated from measurements, 
Zi is a knot of a regular grid on a reference sphere 2p 
of radius R, r is a geocentric position vector, |r| > R, 
j is the resolution level, which is defined by the spatial 
resolution of measurements. The number of the scaling 
coefficients N depends on the density of the regular grid 
and on the size of the modeling area (see Fig. 1) as the 
SSFs are spatially localized (decrease with distance from 
r). 

The SSF is harmonic outside the sphere 2 pg and expressed 
via the Legendre polynomials P, as (Freeden and Michel 
2004) 


oTo 


1 oo R n+l 
; (zi, r) = T X On +1) (=) bjn Pazi r) (8) 
n=2 


where r = |r|, z;, F denote the unit vectors, bi is 
the Legendre coefficient, which determines the spectral 
behaviour of the SSF. We use the Abel—Poisson SSF, 
which is defined by bj, = b”, b; = e7'/?’. The Abel- 
Poisson SSF has a closed-form expression (Freeden and 
Michel 2004), so there is no necessity to calculate the 
sum in Eq.(8). Figure 2 shows the spectral content of 
the Abel—Poisson SSFs of different resolution levels j. A 
greater value of j corresponds to a larger bandwidth of 
the SSF in the spherical harmonic domain and therefore 
to a higher spatial resolution of the gravity model 
Eq. (7). 

The gravity disturbance vector can be represented via 
the scaling coefficients by differentiating Eq. (7) given the 
formula (Torge 2001) 


g(r) = grad T(r) 


as follows 


N 
ég(r) = Yai grad @ ; (z;,1r). 


i=l 


(9) 
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Fig. 1 The ground track of the aircraft’s flight path and the grid of the 
knots of the scaling coefficients 
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Fig. 2 Legendre coefficients b’; of the Abel-Poisson scaling function 
for different resolution levels j 


Assuming that each scaling coefficient a; is constant over the 

time during the aircraft’s flight, the following equations hold: 

a4,=0, i=1,...,N. (10) 

Thus, we obtain the state-space model, which is com- 

pletely described by Eqs.(2)-(6) and Eq.(10) given the 
gravity spatial model Eq. (9). 


The Kalman Filter in the Information Form Let us assume 
that the random errors q „, 4 fr Fas Ef, Cv in Eqs. (4)-(6) are 
zero-mean white noises with known variances. Transforming 
the state-space model Eqs. (2)-(6), (9)—-(10) to the discrete- 
time form, we can further use a Kalman filtering framework. 
The system’s state vector x, at the time moment tk, k = 
0,..., K, where K is the total number of observations, 
includes N + 12 unknown variables 


bvx, Ax, bo, bf, ai, <- AN. 
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Using Kalman filtering, an optimal (in the sense of the 
minimum mean squared error) estimate of the state vector 
is obtained. 


Since the state vector x; at each time moment t includes 
the scaling coefficients for the entire observation area, the 
estimation problem is ill-conditioned. For this reason, we use 
the Kalman filter in the information form (Kailath et al. 2000; 
Park and Kailath 1995). The filter provides the information 
matrix Q; = p and the estimate of the information vector 
up = Pe xy at each time moment t% given the initial values: 
Qo = 0, uy = 0. Here Py denotes the covariance matrix of 
the state vector estimate error. The estimates of the scaling 
coefficients are obtained at the last time moment tx from the 
state vector estimate Xx, which is obtained by solving the 
equation Ox*x = tx. Here, the information vector esti- 
mate ûx and the information matrix Q x are provided by the 
filter at the last time moment tg. As the information matrix 
Q x is ill-conditioned (due to extension of the observation 
area to reduce edge effects as shown in Fig. 1, and due to the 
inverse problem of determining the disturbing potential from 
the gravity disturbance vector), regularization is required. 
Using Tikhonov regularization with a unit matrix J, we 
finally obtain the state vector estimate and the covariance 
matrix of the estimate error, respectively, as 
Px = (Ox +A! (11) 


XK = Prix, 


where A > 0 is the regularization parameter. 


4 Results 


Data Simulation The developed approach was applied to 
processing simulated airborne data. Error-free measurements 
of the inertial sensors and error-free GNSS data (positions 
and velocities) were simulated using the software from Bog- 
danov and Golovan (2017). The simulated aircraft’s flight 
path contained 4 lines flown in the south-north direction 
(Fig. 1) at a constant height of 1,600 m above the reference 
ellipsoid. The initial geodetic latitude is 0.3°. The line 
spacing is about 9.5km, the length of each line is 100km. 
The ground speed at each line is constant and equals 100 m/s. 
The attitude angles are constant at each line (the pitch 
and roll angles equal 0°). The gravity disturbance vector at 
the aircraft’s flight path is synthesized using EGM2008 up 
to d/o 1,500 (equivalent to minimum wavelength of about 
26 km). 

The inertial sensor errors with characteristics shown in 
Table 1 were simulated and added to the error-free data. 
The gyroscope errors were obtained from measurements of 
the calibrated FOG gyroscopes by Optolink recorded during 
a static test at constant temperature. A small constant bias 


V. S. Vyazmin 


Table 1 Statistics for the simulated inertial sensor errors 


Gyrosopes Accelerometers 
Bias instability 0.002—0.004° /h 0.01 mGal/ V Hz 
Noise density 0.1°/h/ V Hz 1.5 mGal/y Hz 


of 0.001°/h was added to the simulated gyroscope measure- 
ments. 

The accelerometer measurement error was simulated as 
a sum of a zero-mean white noise and a bias drift modeled 
as a random walk process (see Table 1). The constant 
bias of 2mGal was added to the horizontal accelerometer 
measurements. The bias of the vertical accelerometer was 
assumed to be determined and removed by comparing pre- 
and post-flight static measurements with a terrestrial gravity 
value. 

The simulated INS measurement frequency is 100 Hz. 
The GNSS velocity error e, is modeled as the time derivative 
of a zero-mean white noise at 20 Hz, the standard deviation 
(STD) of e, is 2-3 cm/s. 


Data Processing The data processing procedure consisted 
of four steps. First, the INS mechanization equations were 
solved using data from the entire flight. The gravity vector 
y, Was associated with the reference gravity field repre- 
sented as the sum of the normal gravity field and a long- 
wavelength part of the anomalous field provided by a global 
gravity model (EGM2008 up to d/o 100 was used). The 
residual gravity disturbance vector (i.e., the gravity distur- 
bance vector with a long-wavelength part subtracted) to be 
estimated varied from —15 to 15mGal for the horizontal 
components and from —20 to 10 mGal for the vertical com- 
ponent. 

Second, the spatial model of the residual disturbing 
potential was designed using Eq. (7). The resolution level 
of the Abel—Poisson SSFs was chosen equal to 8. The 
10km x 4.5 km latitude-longitude regular grid of the scaling 
coefficient knots was used (Fig. 1). The total number of the 
scaling coefficients N equals 154. The dimension of the state 
vector in the Kalman filter is 166. 

Third, the scaling coefficients of the gravity model and 
the INS systematic error parameters were estimated using 
the Kalman filter in the information form and regularization 
Eq. (11) applied at the last time moment tx. The value of the 
regularization parameter A was set equal to sı107!4, where 
sı is the maximum singular value of the information matrix 
Q x, which was found to be suitable. 

Finally, the scaling coefficient estimates obtained at the 
last time moment tx were used to build the along-path 
estimates of the three components of the residual gravity 
disturbance vector using Eq. (9), and the long-wavelength 
part of the gravity disturbance vector was readded to the 
along-path estimates (EGM2008 up to d/o 100). The final 
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along-path estimates of the gravity disturbance vector com- 
ponents are shown in Fig. 3. The estimate errors are shown 
in Fig. 4. 


Numerical Results The statistics for the errors of the along- 
line gravity estimates are given in Table 2. The standard 
deviation varies from 0.6 to 1.9 mGal for the east component 
and from 0.8 to 1.3mGal for the north component. The 
average STDs are 1.2 and 1.0mGal for the east and north 
component, respectively. The mean values are less than 
1.7 mGal for the east and 0.7 mGal for the north component. 
A slightly lower accuracy of the east component, in our 
opinion, is due to the south-north direction of the lines and 
the small total number of the lines. It should be noted here 
that in the case of a smaller number of processed lines, the 
estimation accuracy for the east component is significantly 
worse than for the north component. 


Table 2 Statistics for the errors of the gravity disturbances estimated 
using the spatial modeling of gravity (mGal) 


The accuracy of the along-line estimates of the vertical 
component is better than 0.5 mGal (10), the mean value is 
not greater than 1.5 mGal. 
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Table 3 Statistics for the errors of the gravity disturbances estimated 
with the use of gravity modeling in the time domain (mGal) 


Line no. East North Up 

1 Mean 7.26 0.77 —0.06 
Std 2.64 2.91 0.27 

2 Mean 5.97 5.75 0.76 
Std 2.63 2.11 0.22 

3 Mean —0.35 6.85 —0.41 
Std 4.29 1.60 0.28 

4 Mean —2.85 3.36 1.58 
Std 2.82 4.61 0.31 

Average std 3.09 2.80 0.27 


Comparison with Gravity Modeling in the Time Domain 
The gravity disturbance vector was also estimated using 
the standard approach based on modeling gravity in the 
time domain. Each component of the gravity disturbance 
vector was represented using the Fourier expansion of degree 
15 (equivalent to 25km wavelength). The total number of 
the unknown coefficients (amplitudes of the trigonometric 
functions) equaled 93. The estimation results are shown in 
Fig. 3 and Table 3. The accuracy of the horizontal component 
estimation is significantly lower than in the case of the spatial 
gravity modeling. Large biases and drifts can be observed in 
the gravity estimates. The average STDs are 3.1 and 2.8 mGal 
for the east and north component, respectively. 

For the vertical component, the estimation accuracy (10) 
is 0.3 mGal that is better than in the case of the developed 
approach. The better accuracy is probably due to a better 
approximation of the band-limited gravity signal using the 
Fourier expansion than using the non-bandlimited SSF. The 
mean values of the vertical component estimate errors are 
approximately the same for the two approaches. 


5 Discussion and Conclusions 


An approach based on spatial modeling of gravity using the 
spherical scaling functions was developed for the airborne 
vector gravimetry problem solving. The approach assumes 
the use of a navigation-grade INS and the carrier-phase dif- 
ferential mode of GNSS. From the results of simulated data 
processing, the new approach provides significantly more 
accurate estimates of the gravity horizontal components in 
comparison with the standard approach based on modeling 
gravity in the time domain. The results, however, may be too 
optimistic as simplified models were used for simulating the 
GNSS errors and the inertial sensor errors (only the bias drift 
and noise were simulated). Further, a long-wavelength part of 
the gravity vector was assumed precisely known from a low- 
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resolution global gravity model. Despite this, the developed 
approach, in our opinion, can be put into practice in the near 
future. 

There are, however, a number of issues, the study of which 
is of interest and probably will allow to improve accuracy of 
the gravity estimates. Namely, it should be investigated in 
more detail how the accuracy of gravity estimation depends 
on inaccuracies in the inertial sensor calibration results (e.g., 
on the bias level or a long-term drift in the horizontal 
accelerometers). At the same time, possible ways to improve 
the gravity estimation accuracy should be also studied, e.g., 
using a greater number of parallel survey lines, cross-lines, or 
several overlapping flights. In addition, more careful design 
of the spatial gravity model (e.g., selecting the type of the 
SSF, regularization parameter value, etc.) may also improve 
the resulting estimates. 
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Abstract 


The initial information for the development of high-degree models of the Earth’s grav- 
itational field (EGF) are the results of satellite and ground-based measurements. At the 
same time, satellite measurements carry information on the long-wave structure of the EGF. 
Information on the short-wave structure of the EGF can be obtained only on the basis of 
ground-based measurements. Having organized the determination of deflection of vertical 
(DOV) with a resolution of several kilometers, the local structure of the EGF can be restored 
with the highest possible resolution. This can be done using digital zenith camera systems 
(DZCS). They are automated and allow to determine the components of the DOV at the 
point of placement in real time. The article presents the developed measurement technique 
with a DZCS and the results of its tests at various geographical points in the field. The 
proposed technique, unlike the existing traditional technique, allows to evaluate and take 
into account the calibration coefficients of the DZCS in each series of observations. In 
addition, the new proposed technique does not impose requirements on the accuracy of 
rotation of the telescope around the axis in the horizontal plane and the rigidity of the base 
of the DZCS. The test results of the new technique showed that the standard deviation of 
measurements is about 0.1”—0.3”. 


Keywords 


Calibration coefficients - Deflection of vertical - Digital zenith camera - Observation 
technique 


This method is implemented in the traditional technique used 


1 Introduction in existing DZCS (Albayrak et al. 2019; Hirt et al. 2010; Tian 
et al. 2014; Somieski 2008). 
Components of the DOV can be found in several ways: 2. If the values of the components g+, gy of the gravity vector 


1. If the astronomical ®, A and geodetic B, L coordinates 
of the location are known, the DOV components &, 7 are 
calculated as (Torge 2001): 


£= ©-B, 
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g are known, the components of the DOV are found by the 
following formulas (Brovar 1983): 


E = —arcsin () , 4 = — arcsin (=) ; (2) 


A minus sign means that the DOV in the direction of 
increasing coordinates are considered negative. 
Normalize the gravity vector g in the form: 


n = (A — L) - cosB. (1) 
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where 


A 8x a 


o igh e? 


[gl 


With this in mind, formula (2) will take the form: 


E = —arcsin (x), n = —arcsin(,). (4) 
Thus, the task of determining the values of the DOV com- 
ponents is reduced to calculating the components of the 


normalized vector ĝ. 


2 New Proposed Technique 


The main components of DZCS are a telescope, a CCD 
camera, an inclinometer, a GNSS receiver and auxiliary 
equipment. Measurement data with DZCS in each position 
of the telescope are: a frame of the starry sky, geodetic 
coordinates of the location, exposure time of the frame of 
the starry sky, the current tilt of the telescope from the 
readings of an inclinometer and a temperature value that can 
be determined, for example, from the built-in temperature 
sensor in the inclinometer. 
The unknown parameters of the DZCS are: 


1. Orientation angles ọ, 6, y (Euler angles) between the 
inclinometer coordinate systems (CS) and the CCD cam- 
era CS. They can be determined using rotation matrices 
around the axes (Zharov 2006): 


R= R: (gy): Ry (0) - R: (Y). (5) 

2. The scale factors m, and m, and angle £ between the axes 
of the inclinometer. They can be calculated using a matrix 
of the form: 


myx -sine my» sine 0 
M = 0 m, 0 
0 0 1 


(6) 


3. Temperature coefficients k, and k, of the inclinometer 
axes. They can be determined if the temperature change 
AT during the observation is known: 


A =n?" + ky AT, (7) 

n’, = ny +ky,-AT, (8) 

AT = end — To, (9) 

where kx, ky — temperature coefficients, n7°**, n"°“*— mea- 


sured inclinometer readings without taking into account dis- 
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placement due to temperature; To, Tena — temperature in the 
first and last stationary position of the telescope in a single 
series; n’, n’, — corrected inclinometer readings taking into 
account changes due to temperature. 

4. The components g,, ĝ, of the normalized gravity vector 


£. 

5. During measurements, the DZCS is placed at the mea- 
surement point freely, without orientation to the cardinal 
points. In this regard, it is necessary to determine the 
orientation matrix of the CCD sensor in the local CS 
(topocentric horizontal), in which the OZ axis is aligned 
with the normal to the ellipsoid, the OX axis is directed 
to the north, and the OY axis to the east. Denote this 
matrix A. Matrix A can be calculated on the basis of 
data on the sizes of the CCD sensor, the frame of the 
starry sky, the star catalog, exposure time, geodetic coor- 
dinates, polar motion parameters and the time correc- 
tions received from IERS bulletins (Murzabekov et al. 
2018). 


Based on the foregoing, the parameters of the DZCS are: 
1. measured and calculated: 


(a) CCD sensor orientation matrix in local CS, A; 
(b) measured inclinometer readings on two axes, 
meas meds. 
Ne yor s 


(c) temperature change during a single series, AT. 


2. unknown (10 parameters; the first eight are calibration 
coefficients): 


(a) orientation angles between the inclinometer CS and 
CCD camera, Q, 0, W; 

(b) scale factors and angle between the axes of the incli- 
nometer, My, My, £; 

(c) temperature coefficients of the inclinometer axes, ky, 
ky; 

(d) the components of the normalized gravity vector g, 
§x>Sy- 


Thus, the vector of unknowns of DZCS is as follows: 


X = (9,0, W,mx, my, €, kx, ky, Sx, Sy). 
vea ay 


3 Model of the Proposed Technique 


Multiply the normalized vector g by the matrix A, and then 
by the matrix R. This will allow to obtain the components 
of g in the inclinometer CS. In order take into account the 
inclinometer parameters, multiply by the matrix M and to 
take into account the shifts caused by changes in temperature, 
add term k - AT. Based on this, the model of the new 
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technique can be represented as follows: 


(nme?) = MxRxA,x8+k-AT, 
(nme), = Mx Rx Apx§+k-AT, 


(nme), = Mx Rx Ay x&+k-AT, 
(11) 


where N — number of measurements (number of sta- 
tionary positions of the telescope, N > 10); n”?! = 
(n” od ny öd j” — modeled inclinometer measurements along 
the axes OX and OY; M — matrix for estimating inclinometer 
parameters (6); R — rotation matrix from the CS of the 
CCD sensor to the CS of the inclinometer (5); A — CCD 
sensor orientation matrix in the local CS; $ — normalized 
gravity vector in local CS (3); k = (kx ky j“ — temperature 
coefficients of inclinometer axes; AT — temperature change 
during a single series. To evaluate all unknown model 
parameters, after processing the measurement data in a 


single series, the objective function is formed: 


E-E (EEA E 
eraan), 


where N — number of measurements; (nea) 3 (ay) = 
measured inclinometer readings in the i-th stationary position 
of the telescope, recounted in the projection of the nor- 
malized gravity vector (sines of the respective inclinometer 
reading angles). All unknown parameters are estimated by 
minimizing objective function (12). Nonlinear optimization 
starts with initial values: 


x= (0.0,0.1,1,5,0,0,0,0). 


Least squares optimization is performed using the 
Marquardt-Levenberg method with the numerical calculation 
of derivatives. At the same time, all parameters of the 
model are evaluated simultaneously, i.e. DZCS “auto- 
calibration” occurs in each series. After estimation Êx, Êy, 
DOV components calculated with (4). 

Thus, the new observation technique with the DZCS 
involves obtaining frames of the starry sky, the values of 
the inclinometer in a single series of measurements at dif- 
ferent tilts of the telescope and temperature change dur- 
ing a single series. Measurements in each series can be 
performed in arbitrary directions of the optical axis of the 
telescope and at arbitrary angles in the horizontal plane 
and differ from series to series. The main requirement for 
applying the new technique is the rigidity of the telescope 


- CCD camera - inclinometer system. The accuracy of the 
rotary device and the rigidity of the horizontal plane of 
the DZCS’s location do not affect the accuracy of the final 
values. 
The algorithm of the new technique is presented in Fig. 1. 
The algorithm of the new technique includes three 
stages: 


1. Measurements: obtaining a frame of the starry sky, deter- 
mining the exposure time, determining the current tilt of 
the telescope, measuring the temperature and determining 
the geodetic coordinates in each i-th stationary position of 
the telescope. Moreover, the star catalog, time corrections 
from IERS bulletins and polar motion parameters are 
known in advance. 

2. Data processing: finding stars in the image and deter- 
mining the coordinates of their centers, identifying stars, 
determining transformation parameters and calculating 
the orientation matrix A for all frames of the starry sky. 

3. Calculation of DOV: estimation of the parameters of the 
measurement model in accordance with (12) and calcula- 
tion of the values of the DOV components in accordance 
with (4). 


One the of ways to implement the new proposed technique is 
presented in Fig. 2. 

In this case, the observation process with the new tech- 
nique involves rotating the telescope, CCD camera and 
inclinometer around the vertical axis with a fixed number of 
steps in the horizontal plane twice: with the “initial” zenith 
angle ¢; (first observation cycle) and with set zenith angle 
€2 (second observation cycle). The “initial” zenith angle is 
understood as the state of the initial alignment of the DZCS 
in the horizontal plane according to the readings of the 
inclinometer. 

The advantages of the proposed technique, compared with 
the existing traditional technique, are: 


1. In each series of measurements, a simultaneous assess- 
ment and accounting of all the calibration coefficients of 
the DZCS is made, i.e. there is an “auto-calibration” of the 
device. This avoids additional errors caused by changes in 
calibration coefficients between series of measurements. 

2. Auto calibration process improves measurement 
efficiency. 

3. There are no requirements for the accuracy of rotation 
angles when rotating the device around an axis in the 
horizontal plane. This simplifies the design of the DCZS. 

4. The requirements for ensuring the rigidity of the base 
on which the telescope is located have been reduced. 
It is necessary to exclude any impact on the device 
during measurements in a single stationary position of the 
telescope (approximately 8—10 s). This allows to make 
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Stage 2. Data processing 


Stage 3. 


Determining the coordinates of the 
centers of stars 


Calculation of DOV 


Assessment of 
unknown 


Geodetic 
coordinates 


the star catalog 


parameters of the 
measurement 
model 
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determination of transformation 
parameters 
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Calculation of the orientation matrix 
Fie A of the CCD sensor in the local CS 
x t 


readings 


meas 
3 4 


z3 


Fig. 1 The algorithm of the new observation technique 


Telescope mount rotation axis 


First measurement cycle 
i-th stationary 
position of the telescope 


Second measurement cycle 


Fig.2 An example of the implementation of a new measurement 
technique 


observations on any solid foundation (dirt, asphalt roads 
and sites). This is especially important when measuring in 
the field. 


4 Research of a New Technique 
4.1 Research of Change of the Calibration 
Coefficients 


Tests of the new measurement technique in the field were 
performed on a DZCS (Murzabekov 2017) at five various 
geographical points during 16 stellar observation nights. The 
average number of stars observed in the frame is about 100. 
At each point, at least six series of measurements were 
performed. According to the research results, it was found 
that the number of stationary positions of the telescope 
equal to 24 (12 positions for each rotation) is sufficient. 
Their further increase does not lead to an increase in accu- 
racy. 

In the traditional measurement technique, calibration 
coefficients determined before the start of measurements are 
used as constant values for observations. Studies have been 
conducted to evaluate changes in calibration coefficients 
between series and the effect of these changes on current 
values of DOV. 

An example of the values of the calibration coefficients my 
and m, (inclinometer scale factors) for each series and their 
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Fig. 3 Values and changes of = 1.020 
calibration coefficients m, and m, 
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for each series during * 
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change between the series during the observation at one point 
are presented in Fig. 3. 

As can be seen from Fig. 3, there are a changes in the 
coefficients m, and m,. Consider how the change in the 
coefficients m, and m, between the series affects the DOV 
values for each series and the average values for all series. 
Figure 4 shows the changes in the calculated values of DOV 
depending on the series number (curves | and 2) and the 
change in calibration coefficients my and m, (curves 3 and 4). 

As can be seen from Fig. 4, there is a clear correlation 
between A£ and An and the change in the coefficients m, 
and my. In this case, the differences of single series can reach 
up to 0.4”, the average AE = 0.10”, and An = -0.01”. 

Thus, an uncontrolled change in calibration coefficients 
leads to a shift in the values of DOV, i.e. to the 
appearance of an additional calculation error in DOV. This 
confirms the need to clarify the values of the calibration 
coefficients in each series during observations at each 
point. 


Research of the Influence of the Choice 
of Methods for Processing 
Observational Data on the Accuracy 

of DOV 


4.2 


In the process of research of a new technique we reviewed: 
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e three most used methods for determining the coordi- 
nates of the centers of stars: point spread function (PSF), 
method for approximation of the shape of a star with a 
paraboloid (MAP) and method of segment center determi- 
nation (SCD); 

e four high precision star catalogs: Tycho-2, UCAC4, 
PPMXL, GAIA-DR2; 

e four transformation methods: affine and polynomial 2nd, 
3rd and 4th degree. 

The impact of the choice of processing method on the 

accuracy of the DOV is estimated. The total impact does not 

exceed 0.03”. 


4.3 Measurement Model with an DZCS 
A measurement model with an DZCS can be represented as 
follows: 


h (£, B, ®,n™, xp, yp) = 0, 


| h (nL. A,n™*°, 1y7c,Xp, YP) = 0, OS) 
where tyrc — exposure time of the frame of the starry 
sky in the Universal Coordination Time (UTC); xp, yp — 
current polar motion parameters. In accordance with the 
measurement model (13), the formulas for calculating the 
measurement error of the DOV components with an DZCS 
are written in the following form: 


(14) 


p? 
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Fig. 4 Difference & and ņ, calculated with and without evaluation of calibration coefficients and change of my and my 


where K = 1.1 (with a confidence level of P = 0.95); 
ci — sensitivity coefficients for each component of the error 
(i= 1...11); m; —j-th component of the error (j = 1...9). 

In more detail, the values of each error, sensitivity coeffi- 
cients and calculation examples for several points are given 
in work (Murzabekov et al. 2021). 

For example, for a point on the territory of FSUE VNI- 
IFTRI, the errors of the DOV components according to 
formulas (14) are: mg = 0.36”, m, = 0.24”. Differences in 
errors are due to the dependence of the sensitivity coefficients 


for the DOV component in longitude on the latitude of the 


measurement point. 


4.4 Test Results 


The test results of the new technique are presented in Fig. 5. 
As can be seen from Fig. 5, the value of the standard 
deviation for determining components of DOV is in the range 


0.1”-0.3”. 
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Fig. 5 Standard deviation for determining DOV at five geographical points during 16 observing starry nights 


5 Summary and Conclusions 


Thus, a new technique for performing observations 
with DZSC was developed. It provides the ability to 
“auto-calibrate” the parameters of the device during the 
measurement session in each series. In addition, the 
proposed measurement technique does not require the 
installation of a special hard measuring concrete base 
and precise rotation of the telescope around a vertical 
axis. 

According to the results of testing the device in field 
conditions, the standard deviation was obtained in the range 
of 0.1”—0.3”, which is at the level of world analogues. 
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Abstract 


Long-term recordings of temporal gravity variations were made with the gPhoneX # 
117 spring gravimeter at the Center of Geodesy, Cartography and SDI (TsNIIGAik). 
The purpose of the observations is to obtain reliable parameters of the Earth tides and 
non-tidal component for a site located within the territory of Moscow. At a later stage, 
these parameters will be used to process absolute gravity observations. One year of 
recordings were chosen (April 1, 2016—March 03, 2017). Preliminary data processing is 
fully automated and runs in a program written in Python. In the course of the research, 
the parameters of the instrument drift were determined by the method of piecewise linear 
approximation of the observations. The estimated parameters of the Earth tides were 
obtained according to the algorithms of the ETERNA 3.4 program (Wenzel, Bulletin 
d’ Informations des Marées Terrestres 124:9425—9439, 1996). The non-tidal gravity changes 
(the difference between the experimental and local model data) were compared with the 
model of temporal variations of the gravity field caused by atmospheric loading. Theoretical 
values were provided by EOST Loading Service. 

The studies presented here, are unique for Moscow site and they are one of the first 
attempts to use high-precision recordings of spring gravimeters while computing observed 
tidal and non-tidal parameters. 


Keywords 


Atmosphere loading estimates - Earth tides - Moscow - Non-tidal 


gravity variations 


Gravity variations - 


Earth tides, ocean and atmospheric tides, pole tides, changes 
in groundwater levels, non-tidal loading deformations, etc. 
An algorithm for identifying, accounting, and forecasting 


1 Introduction 


High-precision observations of gravity play a principle role 
in the modern study of the figure of the Earth and its external 
gravity field. It is important to know not only the absolute 
value of gravity at gravity network stations, but also the 
nature of the gravity variations over time. Due to the complex 
interior structure of our planet, gravity varies according to a 
complicated rule. Variations include phenomena such as the 
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all of these factors is developed in each case separately. There 
are some studies for the Earth tides modeling taking place 
in Russia (Spiridonov et al. 2015) but there are very few of 
them based on gPhone gravimeter observations (Antonov et 
al. 2016). In order to obtain and predict the local parameters 
of the Earth tides and analyze loading estimates at the TsNI- 
IGAIK site (ọ = 55.8550°, 4 = 37.5160°), regular observa- 
tions of temporary gravity variations have been made using 
a modern gPhone gravimeter. Such studies have never been 
taken place at Moscow site. It is very important to research 
gravity variation parameters at TsNIIGAiK site thus it is the 
main Fundamental Site of Russian State Gravity Network. 
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2 Data Pre-Processing 


The gPhoneX # 117 spring gravimeter is installed on a pillar 
on the ground floor of the Center for Geodesy, Cartography 
and SDI (Russia, Moscow, Onegskaya St., 26). Continuous 
gravity measurements are being taken from November 2014 
to the present. The interruptions of the recordings have been 
taken place due to power outages and calibration of the 
gravimeter. For this experiment, a continuous 350-days data 
series was selected (from April 01, 2016 to March 03, 2017). 
Frequency of measurements is 5 Hz. The output data used 
for the analysis are: averaged values of gravity variations, 
ambient temperature and pressure, levels positions with 1 s 
resolution. 

Data pre-processing was automated in the Python 3 lan- 
guage and it consisted of eliminating spikes, noise and 
changing the sampling rate to 1 measurement per hour. 
Spikes with deviation greater than 5 jtGal were replaced 
by linearly interpolated values. Noises were cut by apply- 
ing the Savitzky Golay filter (Savitzky and Golay 1964). 
Discretization was made taking into account the Nyquist- 
Shannon sampling theorem (Kotelnikov 1933). 

The polar motion effect on gravity was calculated by 
the Formula, specified in the IAGBN Absolute Observa- 
tions Data Processing Standards (1992). The Earth rotation 
parameters are freely available on the website of the Interna- 
tional Earth Rotation Service (ftp://hpier.obspm.fr/iers/eop/ 
eopc04/). 

The gravity sensor of the gPhone system consists of a low 
drift, metal, Zero Length Spring system (Microg LaCoste 
2012), whose characteristics are determined by the mechan- 
ics of elastic material. Thus the gPhone spring gravimeter is 
characterized by a stable quasi-linear drift of small magni- 
tude (in comparison with quartz gravimeters). The process 
of fitting instrumental drift is made as described by Kang et 
al. (2011). First, according to the characteristics of varying 
trends of gravity residuals time-series (gravity residuals after 
corrections of the solid Earth tide, ocean tidal loading and 
pole tide), the entire time-series are set into segments, then 
segmental polynomial fitting is applied to each segment, and 
finally combined together to compensate entire period for a 
long term gravimeter drift (Kang et al. 2011). 


3 Tidal Analysis 


Let’s compare the differences in the experimental drift cor- 
rected data series with several models of the Earth tides. 
These differences will be called hereinafter a local gravity 
variation component. Theoretical variations were obtained 
for three models of the Earth tides: the harmonic expansion 
of Tamura (1987) (hereinafter—Tamura) for the solid Earth 
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(recommended for data processing by the manufacturer of 
the gravimeter (Microg LaCoste 2012)), tide model WDD 
(Dehant et al. 1999) for an elastic Earth model (obtained 
in the program for processing force variations gravity TSoft 
(Van Camp and Vauterin 2005), as well as the model calcu- 
lated in the Russian program ATLANTIDA3.1_2014 (Spiri- 
donov et al. 2015) using models of the Earth’s structure 
IASP91 (Kennett and Engdahl 1991) and oceans FES2012 
(Carrere et al. 2012) (hereinafter—Atlantida). 

On average over the year, the standard and absolute 
deviations of the local gravity variation component are 22.0 
and 122.9 Gal for the WDD model, 14.8 and 82.0 Gal for 
the Tamura model, 7.6 Gal and 43.7 Gal for the Atlantida 
model, respectively. Thus, the Atlantida model describes the 
tidal gravity variations at the TsNIIGAiK site in the best way. 
It is mainly connected with the fact that Atlantida model 
includes the ocean tides loading effect. 

For analysis, we divide the annual interval into 14-days 
periods. During the year, the local gravity variation com- 
ponent has the same character, but different absolute and 
standard deviations. As an example, Fig. | shows a plot of 
the local gravity variation component for the three models in 
one of the two-week periods. 

Absolute and standard deviations show different values in 
different seasons. Table | shows that the largest deviations 
are observed in the summer and winter months. 

Thus, theoretical model of the Earth tides Atlantida can 
be used for some issues at TSNIIGAiK site. However, local 
parameters of the Earth tides are necessary to conduct high- 
precision studies. The calculations were performed in the 
program ETERNA renewed by Schueller (Schueller 2019). 
The parameters were obtained as corrections to the Tamura 
model. Characteristics of the local parameters are a topic 
for another article. Further we will discuss the differences 
between the experimental data and the local model. These 
differences will be called non-tidal gravity variations. 

Here is the comparison of the local gravity variation 
component and non-tidal component. As an example, Fig. 
2 shows a plot in one of the two-week periods. 

It is possible to conclude that the non-tidal component 
does not exceed 2 Gal in the standard deviation and 
15 Gal in the absolute deviations. The use of the local 
model of the Earth tides reduced unaccounted gravity vari- 
ations by 3 times in comparison with using the theoretical 
model of Atlantida. 


4 Atmospheric Loading Effects 
Let’s now consider the non-tidal gravity variations. The 


main factor that makes them up is the so-called atmospheric 
loading estimates. 
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Fig. 1 Local gravity variation component: differences between observed drift corrected gravity variations and calculated gravity variations using 


the Earth tides models Tamura, WDD, Atlantida 


Table 1 Absolute and standard deviations of the local gravity variation component 


periods deviations, Gal deviations, Gal periods deviations, „Gal deviations, Gal 
TOLO I0 25.79 
2(15.04=28.04) 7.64 32.82 
3(29.04= 12.05) 27.96 
4(13.05—26.05) 39.15 
5(27.0509.06) 33.12 
6(10.06=23.06) 4115 
14.06-07.07) 33.22 
8(08.07-21.07) 39.83 
9(22.07—04.08) 33.66 
10(05.08—18.08) 36.37 
11(19.08—01.09) 29.48 
12(02.09=15.09) 26.30 


Temporal variations in the distribution of atmospheric 
density cause changes in the gravitational attraction of air 
masses during local observations. Also, the load on the atmo- 
spheric masses deforms the Earth’s crust and the sea surface. 
It is known that local measurements of gravity are currently 
corrected for atmospheric pressure variations using a stan- 
dard correlation coefficient of —0.3 microGal/mbar, which 
corresponds to the Resolution IAG No. 9 of 1983. However, 
the current accuracy of relative gravity measurements makes 
it possible to identify the real dependences of temporary 
gravity variations not only on barometric pressure, but also 
on the state of the atmosphere as a whole (temperature, wind, 
humidity, etc.). For this purpose, the University of Strasbourg 


has founded the EOST Loading Service, where one can 
find models of the Earth’s gravity field variations, tilts and 
displacements of the Earth’s surface caused by atmospheric, 
hydrological and induced ocean non-tidal loads. There are 
stations that constantly monitor and analyze the influence of 
the state of the atmosphere on gravity variations to maintain 
this service. 

The EOST Loading Service provides simulated time 
series of gravity variations caused by non-tidal atmospheric 
loading effects according to the data on the state of 
the atmosphere in the Moscow region [http://loading. 
u-strasbg.fr/GGP/atmos]. We have chosen the ERAS 
model (Copernicus Climate Change Service 2017) for 
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Fig. 2 Local gravity variation component and non-tidal component 
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Fig. 3 Unaccounted gravity variations after correcting measurements with drift, Polar motion, local tide and atmosphere loading effects (ERAS 


model) 


processing as recommended in some studies (Olauson 
2018). 

Figure 3 shows non-tidal gravity variations at the TSNI- 
IGAiK point and gravity variations caused by atmospheric 
loading deformations computed using the ERAS model. We 


can see that the ERA5 model closely approximates the 
non-tidal gravity variations. Thus, the unaccounted gravity 
variations are reduced to half of the non-tidal variations after 


discounting the atmospheric loading effect, as it is shown in 
Table 2. 
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Table 2 Absolute and standard deviations before and after using 
ERAS model 


Standard deviations, Absolute deviations, 


Gal Gal 

Before After Before After 
April 1.91 1.12 9.36 6.04 
May 1.75 0.93 6.03 4.63 
June 1.56 1.02 6.33 4.86 
July 2.16 1.26 8.77 5.80 
August 2.90 1.28 6.79 6.46 
September 2.16 ‘| 1.76 7.96 8.21 
October 4.14 1.78 12.95 8.24 
November 3.46 1.86 13.30 9.08 
December 3.23 2.22 13.22 11.01 
January 4.41 2.32 17.22 10.50 
February 5.93 2.44 21.30 11.06 
March 2.93 1.49 11.89 8.76 
Average 3.05 1.62 11.26 7.89 
5 Conclusions 


In this paper, we have shown that the Atlantida model is best 
suited for modeling tidal gravity variations at the site located 
in Moscow. However, when researching high-precision grav- 
ity issues, theoretical modeling is not recommended to be 
used. Using local tidal parameters delivered from the analysis 
of l-year gravity data allowed us to improve the results up 
to 3 times. The non-tidal signal component was less than 
15 Gal peak to peak. 

Finally, in gravity monitoring studies, it is recommended 
to apply atmospheric loading corrections, like that provided 
by the EOST Loading Service (ERAS model). Doing so 
reduces by 50% the unaccounted variations of the local tidal 
model. 
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Abstract 


At present, the uncertainty in measurements of free fall acceleration (FFA) by the designed 
absolute ballistic gravimeter (ABG), developed in VNIIM, attains microgal values (1 
Gal = 1 cm/s’). At such a high level of accuracy, the effects of the interaction of a test 
body (TB) falling in the vacuum chamber of the ABG in the geomagnetic field could 
comprise sources of systematic errors. Furthermore, individual units and systems of the 
ABG itself also can be sources of the magnetic field (MF). We report about more rigorous 
calculations of the possible effects than those performed in the past. A consecutive self- 
consistent method of calculation of the desired correction to the FFA has been developed. 
It results in the differential equation, including not only the field itself, but also its vertical 
first and second derivatives and variation of the velocity along the trajectory. The desired 
correction is obtained by its numerical solution on the basis of the parameters of the field 
generated by the magnet of the stepper motor. The derived correction proved to be at the 
level of a tenth part of microgal. 


Keywords 
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measurements. Such effects, being the sources of systematic 
errors, include, in particular, the effects of the interaction of 
a test body (TB) falling in the vacuum chamber of the ABG 


1 Introduction 


At present, the uncertainty in measurements of free fall accel- 
eration (FFA) by absolute ballistic gravimeters (ABG) attains 
microgal values (1 Gal = 1 cm/s?) (e.g., (Vitushkin et al. 
2020)). The general operation principles of laser-interference 
ABG and the principles of absolute measurements of FFA 
with their help are presented, for example, in (Vitushkin 
2014; Absolutnye gravimetry 2017; Niebauer et al. 1995). 
At such a high level of accuracy, it is necessary to take into 
account a number of effects that influence the result of FFA 


V. I. Sheremet was deceased at the time of publication. 


L. F. Vitushkin - D. D. Gidaspov - F. F. Karpeshin (24) - O. A. Orlov - 
I. S. Khasiev 

D.I.Mendeleyev Institute for Metrology (VNIIM), St. Petersburg, 
Russia 


© The Author(s) 2021 


in the geomagnetic field. Furthermore, individual units and 
systems of the ABG itself also can be sources of the magnetic 
field (MF). Those units are the ion (magnetic discharge) 
pump, electric motor used in some ABG designs, and others. 
Estimates of a possible effect of inhomogeneous MF on ABG 
readings have been obtained in the past, e.g. in (Niebauer et 
al. 1995; Murata 1978, 1980) and references therein. 

In (Niebauer et al. 1995), the effect of the interaction 
of a conducting TB with an inhomogeneous MF has been 
estimated in the case of the FG5 type ABG. The effect is 
due to the emerging Foucault currents. The measured fields 
of the ion pump magnet, servo motor and Faraday optical 
isolators happened to be less than the geomagnetic field, 
estimated at 50 uT, and do not create significant effects on 
FFA measurements with this type of ABG. It was proved 
experimentally that there was no any effect of this kind. 
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To this end, MFs of the order of 100 uT were applied to 
the system with a set of Helmholtz coils. This experiment 
showed absence of effects at a level of 1 Gal. It should be 
noted that in FG5-type gravimeters, the free-fall path of the 
TB is about 21 cm and the fall velocity attains 2 m/s. 

In (Murata 1978, 1980), the motion of a TB charged due 
to the supposed contact phenomena in a geomagnetic field 
under the action of the Lorentz force was considered. The 
estimates also show a negligible effect of the geomagnetic 
field on the motion of the TB. As a consequence, gravimeters 
are usually designed without adjusting for the magnetic 
effects of eddy currents. In the present work, we report about 
more rigorous calculations, we have been undertaken for 
the design of the ABG Grot. The correction is calculated 
on the basis of the parameters of the MF generated by 
the magnet of the stepper motor. The derived correction 
proved to be at the level of a tenth part of microgal. To 
a great extent, it can be taken into account by calibrating 
and comparing the gravimeters. However, in some cases, 
the spread can significantly exceed the indicated value. 
This may be the case if magnetized parts are used. 
Therefore, it is important to understand the nature of 
the correction. Without this understanding, the readings 
of gravimeters of different types can lead to incorrect 
results. 

Bearing this in mind, in the first place the principles of the 
ABG operation are reminded in the next section. The method 
of assessment of the correction is described in Sects. 3 and 4. 
Numerical estimations of the integral correction to apply to 
the determined value of FFA are made in Sect. 5. In the last 
Sect. 6, the results are summarized, conclusions are drawn, 
and further prospects are considered. 


2 Operation Principle of the ABG Grot 


Let us evaluate the uncertainty caused by induction currents 
and their interaction with the geomagnetic field, as well as 
with the MF produced by individual parts of the gravimeter. 
The calculation procedure for all types of MF sources is 
the same. The ABG scheme is usual. The TB falls from 
a small height of 15 cm. Its motion is monitored by the 
laser-interference method. The system generates a set of 
instantaneous values of time t; (zi), i = 1, 2, ..., n, into which 
the TB passes the points z;, fixed in height and spaced half- 
wave apart from each other. Up to the vertical gradient y, the 
FFA at the point z can be written as follows: 

g (z) = go + y (z — z0), (1) 
where go is the FFA value at the point zo. Choosing point zo 
as that where the TB is at the initial moment of time tọ = 0, 
and designating the velocity of the TB at this point vg, we 
arrive at the following set of coupled equations, relating the 
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successive times ¢; with the points z; on the trajectory: 


Zi = Zo + Voli + 


goli” Jy (> 3 


The zo and vg values at the initial time tọ = O form a 
set of initial conditions, which unambiguously defines the 
trajectory. 


3 Physical Premises for 
the Uncertainties Arising in the 
Presence of MF 


Let us consider a simple model. A conductive loop drops 
down city of v(z) along the z axis in an inhomogeneous MF 
with induction B, remaining in the horizontal plane. Since 
the MF is inhomogeneous, the magnetic flux through the 
contour changes during the motion (Fig. 1). As a result, an 
electromotive force appears in the contour, which in turn 
causes circular electric current. As is known, the closed 
frame with current has the magnetic moment m: 


a a =) _ ( =) m 
Ra) = a) 


A 


Fig. 1 Scheme of the fall of the conductive loop in an inhomogeneous 
MF with induction B. The magnetic flux passing through the loop is 
conventionally depicted by the number of lines of force. When the loop 
falls, the number of field lines changes. As a consequence, this induces 
an electric current along the loop and the conjugate magnetic moment 
of the circuit m 
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In Eq. (3), S is the area of the loop, R — its electrical 
resistance. Parameter A characterizes the physical properties 
of the TB: its area, electrical conductivity, and the shape 
in general case. In turn, the produced magnetic moment 
interacts with the external MF. If the latter is inhomogeneous, 
then the loop is drawn into the region with the maximum 
gradient of the MF. Therefore, the force acting on it arises, 
which can be expressed as follows (Landau and Lifshitz 
1960): 

F = —V (mB) = —VU pot, (4) 
where Upot is the potential energy of the interaction of the 
induced magnetic moment with the MF. Dividing the force 
Eq. (4) by the mass of the loop mp, one obtains the desired 
correction to FFA: 


Ao OB, 
Kee L (voB, Z|. 
g mo Oz (: ~ dz ) 


4 Equation for the Motion of the TB 


(5) 


The fall of the TB is described by the second Newton’s law: 


2 


a = g(z) + de.m.(2) 


7 (6) 


In Eq. (6), g(z) is the desired FFA (1). dem, (z) = Ag(z) is 
the correction Eq. (5) taking into account the geomagnetic 
and other electromagnetic fields. The correction for the 
geomagnetic field is small. The correction for the MF of the 
ionic pump can be leveled by the optimal orientation of its 
magnet. Let us examine the correction for the MF created by 
the permanent magnet of the stepper motor. 

Performing differentiation in Eq. (5), one can expand Eq. 
(4) to read: 


A 2 
Ag = -4 (athe + vB Ae + v(%) ) 


=F_+ Fyt Fa 


(7) 


As one can see from Eq. (7), the desired correction 
depends not only on the height z, but also on the velocity of 
the TB at a given point, and moreover, on its own acceleration 
in the point, as 


1 dz 

dv/dz = —~—~. 
a gde 

Therefore, a is included in Eq. (6) twice: both on the left 

and on the right. Thus, the force induced by the MF depends 

on the acceleration of the TB at the given moment and at 

the given point in the trajectory, which is itself determined 


by this force. For this reason, Eq. (6) must be solved in 
a self-consistent way, taking into account this interaction 
of acceleration with the force of electromagnetic induction. 
To this end, we collect both terms containing g = ds 
on the left side of the equality. As a result, after simple 
transformations, one obtains the motion equation for the TB 


in the presence of the MF as follows: 


dz 1 
dt AE Boe is 
t 
AZ [ (aB.\2 32B: 
so- [CE BN. 


By comparing Eq. (8) and Eq. (6), we obtain the desired 
correction to FFA: 


Ag(z) = dem(2) = ag 
_ maT 9) 
ko- (E) + B52 ]} - 8t. 


5 Results 


The total weight of the TB, including the triple prism, 
is P = 58.8 g. The calculation shows that the combined 
coefficient A; = 0.022 m* / Ohm. Let the stepper motor be 
in the x, y plane. It is also necessary to take into account the 
design features of the device, in particular, that its geometric 
center is offset from the origin by approximately 1.5 cm in 
the horizontal plane in both x and y directions. The origin 
is specified by the z axis, which is directed down along the 
trajectory of the stepper motor in the point. The end point 
of the trajectory is 17.2 cm below the stepper motor. In the 
last section of approximately 8 cm length, the trajectory is 
measured using a laser interferometer, and several hundred 
points t; (z;) are used for processing. The measurements have 
shown that the MF of the stepper motor is approximately a 
dipole having components of the magnetic moment: 
M, = 0.16 Am, M, =0.017 Am, M, = 0.014 Am’. 
A numerical solution of Eq. (8) with typical initial con- 
ditions was obtained by the Runge-Kutta method. A repre- 
sentative trajectory z(t) is plotted in Fig. 2a. It is close to the 
conventional parabola. The numerical solution also gives us 
values of instant velocity v(t) along with acceleration 2o, 
Once these values are obtained, one can find instant values 
of the electromagnetic induction correction ae, m.(z), which 
is plotted in Fig. 2a. By fitting the shape of the obtained 
trajectory to the experimental points of the trajectory t; 
(zi), one can find parameters vo, go. The implementation 
of this algorithm depends on a particular device, and also 
on external factors such as seismic vibrations of the base, 
frequency modulation of the laser beam of the gravimeter, 
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Fig. 2 A plot of a part of the 10 
calculated trajectory of the falling 
TB z(t) (a), and the related 


correction to FFA (b), solid 8 
curve, the dashed curve being a 
linear fit 
6 
E 
= A 
N 
2 
0 


10Ag, uGal 


etc. It can be implemented in a number of ways. However, 
in order to evaluate the integral correction to the gravity 
acceleration ae. m,(z), one can proceed as follows. 

Let us use the linear approximation of the curve for 
Ae. m{Z), as Shown in Fig. 2b: 


Ag(z) =k + pz=a+b(z—znm). (10) 
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0,00 0,02 0,06 0,08 


z,cm 


(b) 


In Eq. (10), Zm is the point of measuring the gg value. 
A x? fit has resulted in the following values of the param- 
eters: k = —1.48666 «Gal, and p = 0.194504 uGal/cm. 
First, the latter values provide the following value of the 
desired correction in the lowest point of the trajectory: 
a = —0.011 Gal. Bearing in mind that typical values of 
Ae.m(Z) in the trajectory of Fig. 2 amount to tenth parts 
of microgal, one can assume that such a small final value 
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of the correction might be due to a random compensation 
of individual terms. In any case, one can conclude that the 
desired correction is at the level of 0.1 Gal. Second, using 
the p and k values obtained above, one can derive the value 
of b = 0.1945 uGal/cm from Eq. (10) in a similar way. The 
value of b provides a correction to the gradient of FFA y. At 
the same time, the value of y can be measured independently 
by a relative gravimeter. In principle, this makes it possible 
to compare the present theory with the experiment. 


6 Discussion of the Results 


A consecutive self-consistent method of calculation of the 
desired correction to the FFA has been developed. The 
resulting Eq. (8) explicitly takes into account the fact that 
the correction is determined not only by the field itself, but 
also by its vertical gradient, as well as by the variation of the 
velocity along the trajectory. For this purpose, a consecutive 
consideration of all the three terms F;, F2, and F3 in Eq. 
(7) is performed. The correction is obtained by solving Eq. 
(8) numerically for a representative trajectory of the TB. The 
need for self-consistency is due to the fact that the correction 
depends on the acceleration, on which it has a direct effect. 
It should be noted that in the previous papers (Absolutnye 
gravimetry 2017; Niebauer et al. 1995; Murata 1978, 1980) 
only the term F3 was taken into account, and self-consistency 
was not considered. The obtained estimates indicate that 
the desired value of correction is within the tenth part of 
microgal. 


The developed method also makes it possible to calculate 
the correction to the gravity gradient at the geographical 
point of measuring. In principle, this can serve as a test for 
comparing the theory with the experiment, after the gradient 
has been measured using a relative gravimeter. 
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Abstract 


Volcanic activities sometimes involve gravity changes, and this research is intended to 
establish an observation network surrounding an active volcano using compact absolute 
gravimeters. To simplify the configuration of absolute gravimeters, they are preferably 
operated with a light source distributed from a telecom band (wavelength of 1.5 um) 
laser through optical fibers. To evaluate the accuracy of the absolute gravimeter with the 
telecom band laser, we conducted observations using a prototype gravimeter (TAG-1) with 
frequency-stabilized lasers at both 1.5 jum and 633 nm, and compared these results with the 
expected gravity at the site. Initially, both results showed offsets —187 Gal and —9.6 Gal 
for the 1.5-um laser and the 633-nm laser, respectively (1 Gal = 1078 m/s’). By correcting 
the systematic errors of the photo detectors measured by the synthetic chirp signal, the 
obtained absolute gravity was consistent with the expected value for both wavelengths; 
offsets from the expected gravity were reduced to 6.6 Gal and 5.4 Gal for 1.5 jm 
and 633 nm, respectively. We also evaluated the errors associated with long-distance 
transmission of the 1.5-um laser using a reeled optical fiber (26 km) and an optical amplifier 
and found no degradation in the gravity data from the case of short transmission (10 m). 
These results allow networking of compact absolute gravimeters connected by telecom 
optical fibers that are operated using a common laser and expansion to volcanic areas to 
monitor the gravity change associated with volcanic activities. 
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of magma that is monitored for the prediction of volcanic 
1 Introduction eruptions and to evaluate the transition of volcanic activities. 


Both types of gravimeters, relative gravimeters and absolute 


Volcanic activities often involve earthquakes, crustal defor- 
mation, magnetic anomaly, and other geophysical phenom- 
ena. Gravity change is a direct probe of the mass movement 
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gravimeters, have been used, and the latter can measure long- 
term gravity changes without any instrumental drift with ref- 
erence to an accurate wavelength of the frequency-stabilized 
laser. However, owing to the complex mechanism, large size, 
and high cost, absolute gravimeters have not commonly been 
used for volcanic observations. This research is intended 
to establish a network of compact absolute gravimeters for 
volcanic observations. 

To construct this observation network, absolute gravime- 
ters are preferably operated with telecom band (wavelength 
of 1.5 um) lasers distributed to each gravimeter via optical 
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fibers because conventional lasers (wavelength of 633 nm) 
cannot be transmitted to distant sites because of the loss 
of optical fibers; which is typically 15-30 dB/km for 633- 
nm light, and 0.2-1 dB/km for 1.5-um light. To evaluate 
the accuracy of the absolute gravimeter with a telecom band 
laser, we conducted observations using a prototype gravime- 
ter (TAG-1) with frequency-stabilized lasers at both 1.5-um 
and 633-nm wavelengths, and compared these results with 
the expected gravity of the site. 


2 Gravity Change Associated 
with Volcanic Activities 


Sakurajima is one of the most active volcanos in Japan. A 
devastating eruption occurred in 1914, and small eruptions 
still continue. It has been determined that the amount of 
magma in the magma chamber beneath the mountain is 
coming to that of 1914. Okubo et al. (2013) observed gravity 
changes associated with the volcanic eruption of Sakura- 
jima using an absolute gravimeter. The gravity change was 
10 Gal (1 Gal = 1078 m/s?) at a distance of 2 km from 
the crater; this is only one factor larger than the background 
noise level due to local disturbances such as groundwater 
(Kazama and Okubo 2009). Therefore, observations near the 
crater and networking with a number of gravimeters sur- 
rounding the crater will significantly enhance the detectabil- 
ity of magma motion near the source by averaging the local 
disturbances using a number of sensors. 


3 TAG-1 Gravimeter 


Araya et al. (2014) developed a compact absolute gravimeter, 
TAG-1, and we used it for the evaluation of the availability 
of network monitoring of volcanic activities. To evaluate the 
accuracy of the absolute gravity with the telecom band laser, 
we conducted observations using TAG-1 with frequency- 
stabilized lasers at wavelengths of 1.5 um and 633 nm, 
and compared these results with the expected gravity of the 
site. Figure 1 shows a picture and a schematic diagram of 
the TAG-1 gravimeter. It is comprised of a dropper for the 
free-fall mass and a built-in accelerometer for correction of 
seismic vibrations. By applying the built-in accelerometer 
and the small dropper, TAG-1 is compact and transportable 
for observations. 

The laser light is introduced into the optical unit through 
the optical fiber and is incident to vacuum chambers confin- 
ing the free-fall mass (Free-fall mirror) dropper and a ref- 
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erence pendulum (Reference mirror), both of which include 
retro reflective mirrors forming an interferometer. The inter- 
fered light is guided to photo detectors (PDs) through opti- 
cal fibers. TAG-1 uses a quadrature interferometer for the 
displacement measurement of the free-fall mass, and the 
optical phase is calculated from the detected signals (Hey- 
demann 1981; Svitlov and Araya 2014). From the quadratic 
dependence of the displacement with respect to time, the 
absolute gravity can be obtained. Effects of ground vibra- 
tion acceleration on the gravity are corrected using data 
from the build-in accelerometer using the reference pendu- 
lum. 

TAG-1 can be operated at both wavelengths of 1.5 um 
and 633 nm by using the PDs and the optical unit designed 
for each wavelength. InGaAs-type and Si-type PDs are used 
for the wavelengths of 1.5 um and 633 nm, respectively. 
Beam verticality can be adjusted so that the measuring 
laser beam and beam of the auto-collimator reflected on 
a reference alcohol surface are in parallel. For the ver- 
tical adjustment at the invisible 1.5-um wavelength, the 
measuring laser beam is monitored using an IR (infrared) 
viewer. 


4 Observation Using a Conventional 
633-nm Laser and a Telecom Band 
(1.5 pm) Laser 


We performed gravity measurements in a basement room 
of the main building of the Research Institute of Electrical 
Communication (RIEC), Tohoku University, using TAG-1 
operated with both a conventional iodine-stabilized 633-nm 
He-Ne laser (wavelength of 46 = 632.99081163 nm), and a 
telecom band (1.5 wm) laser (^15 = 1,538.803242 nm, Fig. 
2) which was frequency-stabilized using the acetylene linear 
absorption spectrum with a linewidth of 500 MHz (Kasai et 
al. 2016) and whose frequency accuracy is estimated to be 
107° (Nakagawa and Onae 2004); the latter may realize the 
long-distance distribution of the light source and networking 
of gravimeters. The systematic errors in operation for both 
wavelengths were evaluated. 

The free-fall mirror was dropped every 2 min. The 
obtained data were corrected for seismic noise measured 
by the build-in accelerometer. Figure 3 shows the measured 
gravity using the conventional 633-nm laser and the 1.5-j4m 
frequency-stabilized laser. The theoretical gravity variation 
is shown by the red line based on calculated tidal gravity and 
the absolute gravity measured by the relative measurement 
from a gravity reference point, as described in the following 
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Fig. 1 Picture (upper) and 
schematic diagram (lower) of the 
TAG-1 gravimeter. The optical 
unit shown in dashed red line in 
the lower figure can be replaced 
depending on the laser 
wavelengths 
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Fig. 2 The 1.5-j1m frequency-stabilized laser used for the measurements 
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Fig. 3 Gravity (blue dots) measured using the 1.5-um and 633-nm 
lasers 


section. Slow tidal gravity variations were commonly 
observed for both cases, while offsets were apparent. The 
offsets averaged in the periods were —187 \1Gal for the 
1.5-um laser and —9.6 Gal for the 633-nm laser. This 
may be due to the PDs because mechanical and optical 
configurations are essentially the same for both cases, and 
the laser wavelengths are well defined. We evaluated the 
systematic error of TAG-1 caused by the frequency response 
of the PDs. 
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5 Systematic Error Evaluation 
of the PDs 


Because the frequency of the fringe signal is almost propor- 
tional to the velocity of the free-fall mass, the response of 
the PD causes a systematic error in the gravity measurement 
(Niebauer et al. 1995). To evaluate the error directly, synthet- 
ically modulated laser light that simulates the interferometer 
fringe was applied to the PD, and the difference in obtained 
gravity values from the measurement and from calculations 
could be regarded as estimates of the systematic error. The 
free-fall mass in gravity, g, generates a chirp interferometer 
signal with a frequency rate of df/dt = 2 g/d, where i is 
the wavelength of the laser. Therefore, chirp frequency rates 
of 13.07 MHz/s and 30.96 MHz/s produce g = 9.8 m/s? for 
à = 1.5 um, and 633 nm, respectively. The laser intensity 
was modulated using an electro-optic amplitude modulator 
for the 1.5-jtm evaluation, while a laser diode was used for 
the 633 nm light source. In each case, the chirp signal was 
applied using a function generator whose clock and data 
sampling clock were both locked to the same Rb time base. 
For the measurement of the 1.5-um laser, the chirp frequency 
was set to change from 0.1 MHz to 2.6 MHz in 0.2 s, and it 
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Fig. 4 The observed gravity (blue/light blue dots, left axis) and the 
estimated error of the PDs obtained from the synthetic chirp signal 
(green dots, right axis). Observed 1 and observed 2 are calculated from 
the 1.5-um datasets within 0:30-1:30 and 1:30-2:30, respectively, on 
19 December, 2017, as shown in Fig. 3 


gives 12.5 MHz/s and gı5 = 9.6175202625 m/s? for A15. For 
the measurement of the 633-nm laser, the chirp frequency 
was set to change from 0.1 MHz to 6.3 MHz in 0.2 s, and it 
gives 31 MHz/s and go = 9.8113575803 m/s? for Ag. 

In the data processing of TAG-1, the measured displace- 
ment of the free-fall mirror from fp (time from the release) to 
to + At was fitted using a quadratic function of time, and the 
gravity acceleration, g(fo), was obtained from the quadratic 
coefficient. For small fo, just after a short time from release, 
g(to) showed disagreement with the theoretical dependence 
on fo because the release of the free-fall mirror induced a 
slight recoil vibration; although the total free-fall time was 
approximately 180 ms, the analysis interval was fixed at 
At = 80 ms, and g(fo) was calculated with changing to to 
assess the effect of the recoil vibration. The same calculation 
was applied to the detection of the synthetic chirp signal 
and the systematic errors of the PDs were estimated by the 
difference in g(t) from the calculated gravity (gis or ge) 
obtained from the frequency rate of the chirp signal. Figure 
4 shows the observed g(fo) (blue and light blue dots, left 
axis) and the estimated error of the PDs obtained from the 
synthetic chirp signal (green dots, right axis). The observed 
data for tọ > 60 ms agreed with the error estimation of the 
PDs. The estimation shows small systematic errors for whole 
to, and smaller fọ gives smaller gravity acceleration; the 
decrease of observed gravity at 1.5 um in Fig. 3, calculated 
with fọ = 10 ms and At = 150 ms, is consistent with this 
estimates of errors of the PDs. 
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Fig. 5 The observed data corrected using the estimated errors for the 
633-nm and 1.5-j1m lasers. Theoretical tides were removed and then 
averaged. The expected level based on the relative measurements were 
referenced to AOB-B, as shown by the red dashed line. To see the 
consistency of the corrected gravity, two data sets for each laser were 
calculated and showed similar results 


The observed data, calculated with fọ > 60 ms and 
At = 80 ms, were corrected using the estimated errors, as 
shown in Fig. 5. Theoretical tides were removed from the 
data and then averaged. In contrast to Fig. 3, observed data 
with the 633-nm laser and 1.5-um laser agreed well after the 
correction; moreover, they were consistent with the expected 
level (red dashed line) based on the relative measurements 
using a LaCoste gravimeter referenced to the Aobayama 
gravity reference point (AOB-B) where the absolute gravity 
has been determined. The AOB-B is located 2.3 km west- 
ward from RIEC (Fig. 6). The mean offsets shown in Fig. 
5 were 6.6 Gal for the 1.5-um laser and 5.4 Gal for 
the 633-nm laser, which were compared with —187 «Gal 
and —9.6 Gal, respectively, without the correction in Fig. 
3 

We also evaluated the errors associated with long-distance 
transmission of the 1.5-um laser through the optical fiber. 
The laser light at 1.5 jum was introduced through short 
(10 m) or long (26 km) optical fibers. In this experiment, 
we used a reeled optical fiber in the laboratory. As shown 
in Fig. 7, the measured absolute gravity did not change and 
showed no degradation even when the laser was provided 
through a 26-km-long optical fiber and an optical amplifier. 
Nevertheless, to estimate errors in a practical system in 
the field, environmental effects on the optical fibers, such 
as vibration and thermal disturbances, need to be mea- 
sured. 
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Fig. 6 Location of the Aobayama gravity reference point, AOB-B 


6 Conclusions 


The compact absolute gravimeter, TAG-1, was successfully 
operated with both 633-nm and 1.5-\1m lasers. By correcting 
systematic errors of the PDs measured using a synthetic chirp 
signal, the obtained absolute gravity was consistent with the 
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(KatahireCampus) daN 
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expected value for both wavelengths; the systematic error of 
1.5-um PDs was estimated to be as much as —190 «Gal 
without the correction. These results can lead to network- 
ing of compact absolute gravimeters connected via telecom 
optical fibers operated using a common laser and can be 
expanded to volcanic areas to monitor the gravity change 
associated with volcanic activities. 


Evaluation of Systematic Errors in the Compact Absolute Gravimeter TAG-1 for Network Monitoring of Volcanic Activities 87 


9.80099 r 


9.80098 


theoretical tides 


Gravity acceleration (m/s) 


9.80097 i 
7 


15 8 8.25 


9.80099 r r 


9.80098 


Gravity acceleration (m/s?) 


9.80097 9 9.25 


Day 


Fig. 7 Measured gravity (blue dots) with the 1.5-um laser through a 
short optical fiber (10 m, upper) and a long optical fiber (26 km) with 
an optical amplifier (lower) 
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Abstract 


For the measurement of the acceleration due to gravity, INRiM developed a transportable 
ballistic rise-and-fall absolute gravimeter, the IMGC-02. It uses laser interferometry to 
measure the symmetrical free rising and falling motion of a test mass in the gravity field. The 
launch system is composed of a moveable carriage fixed to two pairs of springs loaded by 
an electric stepper motor, which vertically throw up a corner cube retroreflector in vacuum. 
The interferometer system is a modified Mach-Zehnder interferometer where the launched 
corner cube acts as the reflector in one of the optical arms of the interferometer and the other 
retroreflector acts as inertial reference during the measurement. However, both systems 
entail some practical problems and uncertainty contributions that need to be reduced. In 
particular, the current launch system might cause beam shear and rotational effects due to 
unavoidable small different loadings of the springs, while the current interferometer system 
poses problems in the alignment of the mirrors, which is a highly time-consuming procedure 
and has to be performed before and, sometimes, during the measurement session. For this 
reason, a new launch system consisting of an electric linear motor which produces a linear 
force along its length, and a modified Jamin interferometer system entailing a simpler 
alignment and a better stability in time, have been designed. This works deals with the 
description of these new systems. 
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41 resolution! and, as such, is listed in the BIPM KCDB.” 
The measurement of g is performed using the reconstructed 
rise and fall motion of a corner cube retroreflector, which 
moves vertically in vacuum. An interferometer system is 


1 Introduction 


Absolute measurements of the acceleration due to gravity, g, 


are performed by absolute gravimeters, traceable to the units 
of length and time through their primary standards. For this 
purpose, INRiM developed a transportable ballistic rise-and- 
fall absolute gravimeter, the IMGC-02, which is recognized 
as the Italian primary standard because of the WGG/04- 
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implemented in order to obtain time and space coordinates 
of the trajectory using a visible laser beam. The interfer- 
ometer measures the distance between a free-falling corner 


‘Resolution WGG/04-41 of the first joint meeting of the CCM WGG 
and SGCAG (26-27 May 2004, BIPM): the first joint meeting of the 
CCM WGG and SGCAG recognized the absolute ballistic method of 
measurement of the acceleration due to gravity as a primary method. 
*https://www.bipm.org/kcdb/cme/search?domain=PHY SICS&areald= 
4&keywords=italy &specificPart.branch=19&specificPart.service= 

41 &specificPart.subService=124&specificPart.individualService= 
404&_countries=1&publicDateFrom=&publicDateTo=&unit= 
&minValue=&max Value=&minUncertainty=&maxUncertainty=. 
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Fig. 1 Schematic drawing and picture of the IMGC-02 absolute gravimeter 


cube retroreflector and a second retroreflector mounted on 
the quasi-inertial mass of a vibration isolation system. The 
gravimeter is composed of five main parts: (1) a launch 
system in a vacuum cylinder, (2) an interferometer system 
connected to an anti-vibrating system, (3) a laser body, 
(4) a photodetector and (5) a supporting frame, as shown 
in Fig. 1. A detailed description of the gravimeter can be 
found in D’ Agostino (2006) and D’ Agostino et al. (2008). 
However, the current launch system, which is composed 
of a moveable carriage supported by two pairs of springs, 
and the interferometer system, which is a modified Mach- 
Zehnder interferometer, introduce uncertainty contributions 
and practical problems that need to be overcome. For this 
reason, new launch and interferometer systems have been 
designed. This work deals with the description of the current 
systems and their future improvements. 


2 The Launch System 


The current launch system is composed of the vacuum 
chamber, the test mass and the launching pad. The launch 
chamber can be approximated to a cylinder with a diameter 
of 14 cm and a height of 65 cm. The base of the vacuum 
chamber is made of stainless steel and is supported by three 
legs equipped with levelling screws, which allow the vertical 
alignment of the chamber. The main part of the vacuum 
chamber is a flanged glass pipe, whose bottom part is fitted to 
the base, whereas the upper part is sealed with an aluminium 


cover. A BK7 glass window (5 cm diameter, 1.27 cm thick, 
/10 flatness, parallelism better than 2 arcsec) is positioned 
at the centre of this cover and allows the laser beam to reach 
the test mass. Connections are sealed by O-rings and, inside 
the glass pipe, a Faraday cage shields the test mass from 
electrostatic charges. The base of the chamber has two arm 
pipes: the first one is connected to a low noise turbo pump, 
the second to an ionisation vacuum gauge. A rotary pump 
carries out the coarse evacuation. A pressure of 1 x 107° Pa 
is reached after about 5 h. 

The launching pad (Fig. 2) can be approximated to a 
5 x 5 x 20 cm? parallelepiped and is tightly connected to 
the base of the vacuum chamber. It supports, throws up and 
catches the test mass at the end of its trajectory. The test 
mass is a corner cube retroreflector, which has the property 
to reflect an incident ray parallel to itself regardless of the 
angular orientation. The corner cube has a mass of 0.15 kg 
and lies horizontally on three supports, which are the vertices 
of an equilateral triangle. These supports constitute the upper 
ending part of the catcher, which is fixed to a moveable 
carriage (0.15 kg) on a vertical linear bearing rail, screwed 
to an aluminium rectangular support. The carriage, in turn, 
is fixed to two pairs of series springs and is retained by 
an electromagnet. The starting horizontal position of the 
test mass, resting on the catcher, corresponds to the best 
alignment of the interferometer, i.e. the best overlapping 
between the test and reference beams. The thrust of the 
system (around 27 N) is given by the two pairs of series 
springs (elastic modulus approximately 1.1 N mm! each), 
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Fig. 2 Schematic drawing (left) and picture (right) of the current launching pad 


which are arranged in parallel and loaded through a screw 
gear driven by an electric stepper motor. The vertical motion 
of the test mass is triggered by cutting off the exciting current 
of the magnet. The impulse is transmitted through the three 
supporting points to the test mass. When the mechanical 
action ends, the body hovers and its trajectory is tracked. The 
resulting total stroke of the mechanical system, before the 
release of the test mass, is around 25 mm. When released, 
the test mass has a speed of about 2 m s™! and travels in 
the vacuum chamber a distance of about 20 cm. The launch 
system is automatically reloaded after the pad catches the 
free-falling test mass. 

Ideally, the corner cube retroreflector realizes a point in 
space. Unfortunately, the corner cube’s centre of mass is not 
coincident with its optical centre. For this reason, the corner 
cube is fitted in an aluminium frame, designed in such a way 
as to move the centre of mass of the assembly as close as 
possible to the optical centre. In this way, any rotation of the 
test mass around its centre of mass should affect neither the 
interferometer alignment nor the measured g value. Despite 


the highest accuracy in the realization of the support, a small 
but unavoidable difference between the optical centre and the 
centre of mass is still present. This entails that beam shear 
effects and rotational effects might occur due to movements 
and rotations of the test mass on the horizontal plane, caused 
by (small) different loadings of the springs which generate, 
during the upward launch, lateral forces and moments on 
the moveable carriage. As a consequence, the two interfering 
beams can translate relative to each other; launches affected 
by a fringe visibility reduction above a fixed threshold 
(usually 10-20% of the maximum intensity) are rejected. The 
rotational effect, instead, introduces a centripetal accelera- 
tion, whose vertical component is added to the local gravity 
acceleration g. Experimental tests (D’ Agostino 2006) show 
that the average rotation of the test mass during the trajectory 
is 10 mrad and the associated angular velocity is 25 mrad s~!, 
given a total flying time of around 400 ms. Considering the 
uncertainty in locating the optical centre and the centre of 
mass as a random variable with a uniform distribution within 
+0.1 mm, the expected uncertainty due to the rotation of the 
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Fig. 3 Schematic drawing (left) and picture (right) of the new launching pad 


corner cube is 3.6 Gal and represents one of the largest con- 
tribution to the uncertainty budget (D’ Agostino et al. 2008). 

To overcome these issues, a new launch system has been 
designed (Fig. 3). It can be approximated to a 7 x 7 x 26 cm? 
parallelepiped and consists of an electric linear motor, which 
has its stator and rotor unrolled to produce a linear force 
along its length, to which the moveable carriage and the 
catcher are fixed. As in the current system, the moveable 
carriage slides on a linear bearing rail screwed to an alu- 
minium rectangular support. In this way, once aligned, beam 
shear and rotational effects should be minimized. A 50% 
decrease in the number of rejected launches is foreseen, and 


the rotational effect uncertainty is expected at least to halve. 
The linear motor has a maximum stroke of 780 mm, a peak 
force of 67 N and a continuous maximum force of 14 N. The 
presence of motor flanges enables the easy mounting of the 
linear motor on the aluminium structure, while the clamping 
plate design enables quick assembly and disassembly of the 
linear motors without disassembling the flange. The rate of 
movement of the magnetic field is electronically controlled 
via software, which allows one to track the motion of the 
rotor. In this way, it is also possible to design the downward 
motion profile in order to collect the test mass at the end of 
its falling motion. 
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3 The Interferometer System 


The current interferometer system is a modified Mach- 
Zehnder interferometer (Germak et al. 2002; D’ Agostino 
2006). A scheme is reported in Fig. 4. The light emitted by 
the He-Ne laser (à ~ 633 nm) passes through an optical 
Faraday isolator (FI) which avoids any backreflection into 
the laser, possibly changing its wavelength and thus affecting 
the measurement. After the Faraday isolator, the beam is 
enlarged to 2 mm spot size by a beam expander (BE) and 
proceeds vertically toward a beam splitter (BS). The reflected 
beam proceeds horizontally toward a small corner cube (Scc) 
used to observe the vertical orientation of the beam through 
a telescope (T). The transmitted beam, vertically directed, 
is deviated by a moveable mirror (M1) which can be rotated 
around two axes to adjust the direction of the shifting arm 
of the interferometer to the true vertical. These elements are 
needed to check the verticality of the laser beam. Afterwards 
the beam enters horizontally into the optical prism (OP) 
situated beneath the reference corner cube retroreflector 
(Rr). The OP has been designed in such a way that one half 
is a reflecting mirror and the other half is a beam splitter. The 
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laser 


The modified Mach-Zehnder interferometer 


Fig. 4 Scheme of the current modified Mach-Zehnder interferometer 
(left) compared to the new modified Jamin interferometer (right). 
Abbreviations: OP optical prism, BE beam expander, FI Faraday iso- 
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incident light is divided into two beams by OP. One of the 
beams, which represents the fixed arm of the interferometer, 
proceeds horizontally and goes straight toward the detector 
(D); the second beam travels vertically down to the test 
mass’s corner cube retroreflector (Tr) and forms the shifting 
arm of the interferometer; it is reflected back to Rr and 
strikes the movable mirror (M2) after being reflected by the 
OP. The adjustment of M3 allows the beams to recombine 
at the second beam splitting portion of OP. The interference 
fringe causes cyclic variations in intensity, which are due to 
the change in the difference between the two optical path 
lengths at every half-wavelength displacement of the test 
mass retroreflector Tr. The detector converts the output 
light of the interferometer to an electric signal. The test 
mass trajectory is measured by timing this electric signal. 
Since the recombining beams have to be coaxial in order to 
avoid distortions on the laser interference fringes, the angular 
position of the mirror Mz has to be adjusted every time before 
the measurement and has to be monitored during the mea- 
surement. The alignment is achieved by rotating the mirror 
Mp around its two axes through a piezoelectric tilt actuator. 


Tp 
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lator, BS beam splitter, D detector, T telescope, IP interference pattern 
control, M mirror, Rp reference retroreflector corner-cube, Tp test-mass 
retroreflector corner-cube, Scc small corner-cube 
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Fig. 5 3D-Scheme of the new modified Jamin interferometer on the horizontal plane 
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Fig. 6 Bottom view of the 3D-Scheme of the new modified Jamin interferometer integrated in the IMGC-02 rise-and-fall ballistic gravimeter 


Unfortunately, this operation entails practical problems, 
is highly time consuming and has to be performed before 
and during the measurement session. For this reason, a 
new modified Jamin interferometer (Shamir 1999) has been 
devised (Fig. 4). Such system is similar to the modified 
Mach-Zehnder interferometer except that the two beams 
directly recombine on the OP, thus the movable mirror M2 


is removed. A 3-D detail of the main part of the new system 
is depicted in Fig. 5, while its assembly in the gravimeter is 
shown in Fig. 6. The alignment of the recombined beams is 
possible by just shifting the reference corner cube retroreflec- 
tor Rp along the horizontal plane. The main advantages are 
a simpler and faster alignment of the two beams and a better 
stability in time entailing an expected setting time reduction 
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of around 25% and a further decrease of rejected throws. 
The new scheme has a potential Abbe error because of the 
removal of Mz, which was used to correct the misalignment 
of the beams traveling in the two arms of the interferometer; 
using high optical quality corner cubes and prism (angular 
accuracy within 1 arcsec and a flatness within 4/10) the 
uncertainty contribution due to the Abbe error is minimized. 


4 Conclusions 


For the measurement of the acceleration due to gravity, 
INRiM developed a transportable ballistic rise-and-fall abso- 
lute gravimeter: the IMGC-02. Currently, the launch system 
is composed of a moveable carriage supported by two pairs 
of series springs in parallel, which are loaded by moving 
down the carriage through a screw gear driven by an electric 
stepper motor. Nevertheless, it is likely that lateral forces 
arising during the upward launch of the test mass, caused 
by unavoidable small different loadings of the two springs, 
make the test mass move on the horizontal plane and rotate. 
Therefore, the two interfering beams can translate relative to 
each other, so that beam shear and rotational effects might 
occur. To overcome these issues, a new launch system has 
been designed. It consists of an electric linear motor, which 
has its stator and rotor unrolled to produce a linear force 
along its length, fixed to the moveable carriage. The rate of 
movement of the magnetic field is electronically controlled 


to track the motion of the rotor. In this way, beam shear and 
rotational effects should be minimized. For what concern the 
interferometer system, at present a modified Mach-Zehnder 
interferometer is adopted. Unfortunately, the alignment oper- 
ation entails practical problems, is highly time-consuming 
and has to be performed before and, sometimes, during the 
measurement session. For this reason, a new modified Jamin 
interferometer has been devised. This system is similar to 
the modified Mach-Zehnder interferometer except that the 
two beams directly recombine on the optical prism, thus 
the movable mirror is removed. The main advantages are 
a simpler alignment of the two beams and better stability 
in time. The new interferometer system should guarantee 
measurements that are more robust and a faster setup of the 
gravimeter. 
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Twelve Years of High Frequency Absolute 
Gravity Measurements at the UK’s Space 
Geodesy Facility: Systematic Signals 

and Comparison with SLR Heights 


Victoria Anne Smith, Graham Appleby, Marek Ziebart, and Jose Rodriguez 


Abstract 


Absolute gravity measurements taken on a near-weekly basis at a single location is a rarity. 
Twelve years of data at the UK’s Space Geodesy Facility (SGF) provides evidence to show 
that the application of results from international comparisons of absolute gravimeters should 
be applied to data and are critical to the interpretation of theSGF gravity time series of 
data from 2007 to 2019. Though residual biases in the data are seen. The SGF time series 
comprises near weekly data, with exceptions for manufacturer services and participation 
in international instrument comparisons. Each data set comprises hourly data taken over 
1 day, with between 100 and 200 drops per hour. Environmental modelling indicates that 
the annual groundwater variation at SGF of some 2 m influences the gravity data by 3.1 
uGal, based upon some measured and estimated soil parameters. The soil parameters were 
also used in the calculation of the effect of an additional telescope dome, built above the 
gravity laboratory, and have been shown to be realistic. Sited in close proximity to the 
long-established satellite laser ranging (SLR) system and the global navigation satellite 
systems (GNSS) the absolute gravimetry (AG) measurements provide a complimentary 
geodetic technique, which is non space-based. The SLR-derived height time series provides 
an independent measurement of vertical motion at the site which may be used to assess the 
AG results, which are impacted by ground motion as well as mass changes above and below 
the instruments. 


Keywords 


Absolute gravimeter - Gravimetry - International Terrestrial Reference Frame (ITRF) - 
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1 Introduction 


Characterisation of key geodetic sites for the improvement of 
their products is important for the demands of geodesy and 
the various reference frames which are built upon the data 
from the worldwide geodetic observatories. The addition 
of absolute gravity (AG) measurement capabilities at the 
Space Geodesy Facility (SGF) was prompted by the growing 
demands of the Global Geodetic Observing System (GGOS) 
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and the European Combined Geodetic Network (ECGN). 
An early goal for the observatory was to determine if the 
geodetic products from the long established satellite laser 
ranging (SLR) and GNSS, could be improved through better 
characterisation of any un-modelled or mis-modelled signals 
by a non-space based AG technique. Determination of the 
glacial isostatic rate for the southeast of the UK was an 
additional long term goal. 

In this paper, the results of 12 years of near-weekly 
AG data are presented. The environmental effects and cal- 
culations are briefly discussed. Emphasis is placed on the 
interpretation of the time series as a whole, where instru- 
mentational inter-comparisons, data offsets and bias offsets 
are discussed and analysed. Finally, a comparison is made 
between the AG time series and station heights derived from 
SLR. 


2 Site Information 


The UK’s Space Geodesy Facility is funded by the UK 
Natural Environment Research Council through the British 
Geological Survey. It is located in East Sussex, less than 
5 km from the English Channel. The geodetic techniques 
of SLR, GNSS and absolute gravity measurements form 
the principle operations at the observatory, though numer- 
ous additional sensor data, including sun photometry and 
automated measurements of groundwater depth, atmospheric 
visibility, pressure, temperature and humidity, are recorded. 
The site is compact with each geodetic technique contained 
within 25 m. An additional GNSS antenna is located around 
100 m distant. The gravimetry laboratory is located SE from 
the SLR, three meters below ground level as measured at the 
borehole, in a semi-sunken bunker that has approximately 
1.2 m of soil above it. Unusually for geodetic sites, the 
subsurface is comprised of clay, with no bedrock beneath 
within 30 m. The local soil type is known to be Weald 
clay (British Geological Survey website!). Absolute gravity 
measurements have been taken at the observatory since 
2006 on a near-weekly basis. The standard measurements 
comprise of hourly data sets with 100-200 drops per set 
taken for 25 h once a week. 


3 Hydrology 


The hydrology surrounding the gravity laboratory is slightly 
complex due both to its semi-sunken nature and the slope 
on which the SGF sits (Fig. 1). It has been well docu- 
mented that water movement should be accounted for when 
analysing gravity data (Makinen and Tattari 1990; Harnisch 


"http://mapapps.bgs.ac.uk/geologyofbritain/home.html?. 
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Gravity 
Laboratory 


Fig. 1 Schematic to show the semi-sunken nature of the gravity labo- 
ratory 


and Harnisch 2006) groundwater, soil moisture and rainfall 
all change the mass around gravimeters and therefore, it 
should be understood in order to interpret the data correctly. 
Groundwater data have been recorded at Herstmonceux since 
the mid 1990s by an automatic data logging system in a 
borehole, measurements are recorded as depth below ground 
level. There is currently no capability for soil moisture to be 
recorded. 

However, using estimates made by testing soil samples, 
approximations of the influence of these hydrologically- 
induced signals on the AG data have been made. The density, 
porosity and specific yield were needed to estimate the 
differences between wet and dry soils both above and below 
the gravimeter. The modified Bouguer plate corrections for 
soil moisture and ground water are: 


885m = 20GH py SP = 4.2H8P (1) 


gew = 20GPpydH = 4.2P8H (2) 


Where 685m indicates the effect on gravity given by a change 
in soil moisture, 5g,, indicates the effect on gravity given 
by a change in groundwater depth, G is the gravitational 
constant, pw is the density of water, ôP is the change in the 
water filled pore spaces in the soil and ôH is effective change 
in depth to groundwater. 

The effect on gravity, due to the variation of water content 
in the 1.2 m of soil above the gravimeter, was calculated to 
be less than 1 Gal. The influence on the gravity data due to 
seasonalvariation in the groundwater height, of between10 
and 12 m below the gravity laboratory, has been calculated 
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Fig. 3 Regimes split by event 
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to be 3.1 Gal at maximum. Application of the ground- 
water correction gives an almost imperceptible change to 
the appearance of the AG time series. Full details of these 
calculations can be found in the PhD thesis by Smith (2018).? 


4 Interpretation of the Time Series 


The weekly gravity data as presented in Fig. 2 clearly 
contain some interesting features. However, these features 
become of particular significance when important events 
such as changes to the instrumentation or environment are 


? Available from the British Library or UCL Discovery website. 


Uncorrected Gravity Time Series 2006 - 2019: Events 
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Gravity Time Series 2006 - 2019 
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2014 2016 2018 


highlighted, as shown in Fig. 2. These events appear to cause 
clear offsets in the time series of AG data. If the data are 
split into slices, dictated by each event, then each resulting 
‘regime’ of data, as shown in Fig. 3, can be analysed and the 
differences between them can be studied. 

It is interesting to note that the data in each regime falls 
within a 10 Gal range but are offset in mean value. Also, the 
standard deviation about the mean of each regime is similar: 
1.17, 1.31, 1.39, 1.56 and 2.28 Gal. 

It should also be noted that the precision of the AG 
measurements has been decreasing for several years, 
which is thought to be driven by the environment, as 
comparison data has not indicated any problem with the 
instrument. 
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Fig. 4 Comparison corrected 
time series 
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If the data are to be interpreted as a whole, it is clear 
that the offset changes, coincident with instrumental visits 
to the manufacturers, need to be accounted for. Since the 
measurement of gravity is dependent on instrumentation, as 
well as spatial and temporal variables, no estimation of the 
accuracy of the instruments is possible. To quantify results 
of absolute gravimeters the best solution is to obtain relative 
offsets between instruments during international gravimeter 
comparisons (Francis et al. 2013). Relative offset values 
obtained from comparison results were obtained three times 
during this time series for the FG5 used (FG5-229). The 
comparisons results for 2007 and 2015 gave offset results 
for the SGF instrument of +1.5 + 0.8and 0.08 + 0.8 Gal 
respectively (Francis and van Dam 2010; Pálinkáš et al. 
2017). A basic mini-comparison was carried out in 2012, 
courtesy of O. Francis, at the Walferdange Underground Lab- 
oratory, with another FG5 (FGS-216) which itself was found 
to have an offset of +1.8 + 3.1 jtGal during the international 
comparison of 2011 (Francis et al. 2013). Since FG5-229 
was found to be in very close agreement (0.2 + 1 Gal) 
with FG5-216, the offset of +1.8 has been used for FG5- 
229. Unfortunately, since there is no supporting comparison 
for the early data, before the SIM (system interface module, 
through which the majority of signals are passed and contains 
control electronics) repair in 2006, the data from the first 
regime shown in Figs. 2 and 3 have been discounted from 
further analysis at this time. 

The offset corrections from the comparisons have been 
applied and the resulting time series is shown in Fig. 4. The 
data from 2007 to 2013 have been significantly smoothed 
and now appears to give consistent data. However, there are 
significant discontinuities remaining thereafter. 


Comparison Corrected Gravity Time series 2006 - 2019 


2010 2014 2016 2018 


The offset in the data from 2016 onwards is of known 
origin; a small telescope dome was built above the gravity 
laboratory. An approximate calculation of mass change, 
based upon volume of clay extracted from above the gravity 
laboratory and the addition of concrete and brick, implies 
an offset in the gravity measurement of 3.6 Gal. This is 
in reasonable agreement with the observed mean offset as 
discussed below. However, the data from 2013 are concern- 
ing. They have a large residual offset after the comparison 
results are applied. We have no evidence to support any local 
environment changes that could have accounted for this level 
of change over the 3 months that the instrument was away 
at the manufacturers. Although heavily researched, this bias 
remains of unknown origin. 

Several methods could be employed to determine the 
biases in the 2013—2016 and 2016-2019 data sets. However, 
using the simplest method, calculating the mean values of 
each regime, proves to be an acceptable method in this case, 
since the difference between the mean values of each regime 
compares favourably with both the offsets obtained from 
the 2007 and 2011 comparisons and the estimated effect 
of mass changes due to the additional telescope dome. The 
comparisons in question gave a total offset of 3.5 Gal, while 
the difference between the mean values is 3.6 Gal. The 
calculated change due to the additional dome built above the 
laboratory, based upon soil samples, indicated an expected 
offset of 3.6 Gal, whereas the difference between the mean 
values gives 4.8 „Gal. It should have noted that the dome 
calculations did not account for the additional masses from 
the telescope, telescope mount, pier or dome, which would 
all increase the expected offset. Furthermore, in earlier work, 
the magnitude of the bias in 2013 was taken to be 7.3 Gal 
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Fig. 5 Time Series corrected by 
comparison offsets and mean 
differences 
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Fig. 6 Campaign simulation, 
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obtained by ‘best-fit’, whilst the difference between the mean 
values gives an offset close to this of 7.2 Gal. When these 
offsets are applied, in addition to the offsets provided by 
comparisons of instruments, the results are shown in Fig. 5. 


5 Campaign Simulation 


To emulate low frequency FG5 data a ‘campaign style’ 
simulation was applied to the SGF time series; by taking one 
data point per year from the full SGF time series. Project 
data, comprising 25 h of data, were selected around the same 
epoch each year and plotted in Fig. 6. Interpretation of these 
results could be critical for a country-wide project: it could 


2010 


Comparison Corrected Gravity Time series 2006 - 2019 


2010 2012 2014 2016 2018 


Year 


Simulation of 'campaign' single project measurement 2006 - 2019: Uncorrected 


2012 2014 2016 2018 


be tempting to apply an offset for the 2013 to 2017 data and 
estimate a glacial isostatic adjustment (GIA) rate from the 
residual. Such an interpretation would be critically flawed 
given our understanding of events in the whole SGF time 
series. 

Bias and comparison corrected results for this campaign 
simulation are given in Fig. 7. 


5.1 AG vs SLR Heights 

The SGF is an International Laser Ranging Service (ILRS) 
Analysis Centre (AC) and is responsible for computing 
reference frame solutions that are subsequently made freely 


102 


V. A. Smith et al. 


Fig. 7 Campaign simulation 
data corrected for bias, dome and 
comparisons 
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Fig. 8 AG data, converted to 
heights, plotted in red, with SLR 
derived station height changes 
plotted in green 
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available to the community via ILRS Data Centres.* The 
SGF AC contributed multi-year solutions towards the real- 
isation of the International Terrestrial Reference Frame (e.g., 
ITRF2014, (Altamimi et al. 2016)) using global SLR mea- 
surements to the LAGEOS and Etalon satellites. In addi- 
tion, research carried out by the AC, (Appleby et al. 2016; 
Rodriguez et al. 2019) after the work that contributed to 
ITRF2014, has identified subtle systematic effects in the 
range measurements at each station of the ILRS network 
that impact in particular the derived station heights at the 
few mm to 10 mm level. As a consequence of this work, 
new solutions for station coordinates that take account of 


$https://ilrs.cddis.eosdis.nasa.gov/data_and_products/data_centers/ 
index.html. 
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these small systematics have been computed and may be 
considered close-to bias free. The resulting weekly-averaged 
heights of the Herstmonceux site thus provide a reference 
height series against which to measure the stability and 
potential systematics in the AG gravity observations; the AG 
measurements will be impacted by height changes, as tracked 
by the SLR solutions, as well as mass changes above and 
below the AG instrument, primarily hydrological changes 
as discussed in Sect. 3. Therefore, the data from the two 
techniques are not expected to match precisely. The SLR 
and AG time series (converted to height by the use of a 
multiplication factor of —5 mm/wGal (Teferle 2009)) are 
shown in Fig. 8, where the data has been aligned in the 
vertical axis by matching the first data point of both the AG 
and SLR data. It is very interesting and promising to note 
that the AG series provides a smoother representation of the 
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height variation of the site than is determined by the SLR 
time series. 


6 Conclusion 


The high frequency, near weekly, AG data from SGF show 
clear evidence of the importance of testing the instruments 
after they have been serviced and after changes to the 
instrument have been carried out by the manufacturer. The 
implementation of the comparison results from international 
meetings is seen to be critical, though not the complete story 
for the bias offset seen in 2013. It is alarming that this bias 
coincides with the major instrument change of the upgrade 
of the FG5 to an FG5X. The results indicate that instrument 
comparisons should be carried out as often as is practical for 
all users of FG5s. 

SLR-derived site height variations are used to validate the 
concept of applying the corrections made to the AG time 
series that results from inter-comparisons, bias estimation 
and new dome building. SGF is currently testing an addi- 
tional site in the UK with an aim to be able to offer a mini- 
comparison site for the community. 
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Abstract 


The paper presents research of the evolution of resolution capabilities and approximation 
accuracy of the last decade global geopotential models based on space gravimetric 
measurements of CHAMP, GRACE and GOCE missions using their spectral characteristics. 

The comparison between the model data with point measurements data on gravity 
anomalies and quasigeoid heights for the territory of Novosibirsk region is shown. Based on 
the research results the conclusion was drawn that accuracy characteristics of current global 
models under test built by the results of satellite gravimetry missions do not achieve the 
specified accuracy of 1 cm and 1 mGal on the territory under investigation. The research has 
made it possible to state that at the current technological and methodical level the potential 
has been reached as concerns EGF models resolution and accuracy enhancement. 
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The object of this paper is the study of evolution of 
resolution and approximation accuracy for global models by 
their spectral characteristics. 


1 Introduction 


The study of the resolution and approximation accuracy of 
the global gravity field (EGF) models is a vital task as 
regards fine structure of the gravitational field. Current mod- 
els in the form of standardized coefficients of geopotential 
spherical harmonics make it possible to construct detailed 
high accuracy digital models of gravity fields characteristics 
(Koneshov et al. 2013; Aka Blush Ulfred 2019; Erol et al. 
2020; Abd-Elmotaal et al. 2020). Reliable estimate of EGF 
models accuracy increases opportunity of their application in 
geodynamics, geophysics and navigation. 
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Estimated approximation accuracy criterion for quasi- 
geoid height and gravity anomalies implies the values of 
1 cm and 1 mGal respectively, as specified by project GOCE 
developers. 


1.1 Theoretical Part 
When modelling the Earth’s gravity field, the attraction 
potential expansion in a Fourier series is used by the system 
of spherical harmonics of geocentric coordinates—radius 
vector r, latitude @, and longitude À in the form (Hoffman- 
Wellenhof and Moritz 2005) 


V(o.A,r) = = [! + =). z (Can cos mÀ 


4S nm sin m)Pan( sini o)] 
(1) 
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where fM—geocentric gravitational constant; 

a -—Earth’s equatorial radius; 

Crm H Snm—normalized dimensionless harmonic factors of 
geopotential; 

P nm (sin y)—normalized associated Legendre polynomials. 


Summing-up in formula (1) is performed to infinity, and 
the Earth’s gravity field models are limited by the maximum 
degree No. Series (1), limited by maximum No presents spec- 
tral expansion of the Earth’s gravity field structure by waves 
lengths 360°/2No that corresponds to the spatial angular 
resolution. 


1 fo} 
A= 0 : 
No 


(2) 


Simulation error of the global Earth’s gravity field V(@, 
À, r) by Fourier series Vy (@, A, r) in point P(g, A, r)Ew with 
approximate harmonic coefficients LC wins Sa mt is written as 


Aw = p(V, Vy) + max | ôN |, (3) 
where 
p(V, Vn) = max | V(P) — Vn (P) |; (4) 


d6y—simulation error of field V(P) due to the error of 
harmonic factors {6 Cnm, 5Snm}. 

Geopotential approximation dispersion error V(q, A, r), 
taking into account the error of harmonic factors, is estimated 
by formula (Kanushin et al. 2015) 


N 
D {An} = D-Y (Dn = dy) 


n=0 


(5) 


where D is a dispersion of the initial potential V with 
unlimited spectrum. The degree dispersion of geopotential 
is calculated by formula 


ps) E Caah 


m=0 


(6) 


Formula (7) is used to calculate the degree dispersion of 


errors 8C,,,, and Sn of geopotential harmonic factors V 


nm 


nl (Tar: 


m=0 


(7) 


Dependence of approximation error ¢ of gravity anoma- 
lies and quasigeoid heights due to series (1) limit in the 
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form of 


2 N 
dpat e 
ety Dey ZNL Da di, 
spherical degree dispersions of gravity anomalies and quasi- 
geoid heights expressed by factors ACnm and Sym 


Das, = 
Ds, 


where ACim = Cam — T? „—the difference between 
coefficients of normalized spherical functions of real and 
normal gravity fields; 

T? „—coefficients of normal geopotential are referred to 
ellipsoid WGS- 84. 

The degree dispersion of the errors of gravity anomalies 


and quasigeoid heights coefficients is written as 


i- 
ds, 

Experimental. The study of evolution of the resolution 
and accuracy for 70 models under test for global Earth’s 
gravity field was conducted by the comparison results of 
geopotential harmonic coefficients degree dispersions and 
their errors as well as the computed (by them) degree 
dispersions of gravity anomalies and quasigeoid heights. 

For satellite models under test their resolution and accu- 
racy results according to their spectral characteristics are 
presented in the time range of 2008-2015 (Fig. 1). At 
the present moment the maximum expansion degree of the 
satellite models has achieved No = 300. The given expansion 
degree is limitary for satellite models built by current space 
gravimetry missions. 

The evolution of the spatial resolution of the models 
built by the last decade satellite data shows that the gravity 
anomalies models resolution averages to 100 km, and for 
quasigeoid heights 200 km. 

The results of global models resolution obtained by the 
combined data with error approximation 1 cm and 1 mGal are 
presented for the period of 1996-2019 (Fig. 2). The degree 
dispersions for the combined models under test practically 
coincide to the 240th degree. This testifies to the fact that 
the low-frequency harmonic segment of these geopotential 
models has been thoroughly studied. 

The data on the evolution of resolution capability for the 
gravity anomalies and quasigeoid heights models obtained 
by the combined data shows that the resolution capability 
was changing insignificantly during the last 15 years, i.e. 


n 


D (ACam + Fim)» © 
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Fig. 1 Resolution capabilities for spectral characteristics obtained according to low-level models, when approximation errors of 1 cm and 1 mGal 
are achieved (period 2008-2015) 


about 180 km for quasigeoid and about 70 km for gravity practically coincide over all the frequency range under study. 
anomalies. This means that the same initial data were mostly used for 

The results of spectral estimates of global geopotential constructing these models. Significant improvements for the 
ultrahigh-degree models show that their degree dispersions resolution of gravity field ultrahigh-degree models (Fig. 3) 
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Fig. 2 The results of calculating the resolution capabilities from the spectral characteristics obtained from the combined data, with approximation 


error of 1 cm and 1 mGal (period 1996-2019) 


built in 2008-2019 have not been observed. For the gravity 
anomalies models the resolution is about 45 km and that for 
quasigeoid height—about 170 km. 

Analysis of the obtained research results as regards evo- 
lution of gravity field models resolution capability (for com- 
bined ultrahigh-degree models made in 2008-2019) revealed 
that it has not practically changed. Resolution of satellite 
models (in the period of 2008-2019) changed for the gravity 
anomalies model from 200 km to 100 km; for the quasigeoid 
model it changed insignificantly, i.e. from 230 km to 200 km. 

Analyzing the spectral estimator changes from one geopo- 
tentional model to another a clear view may be formed 
of the resolution capability of certain harmonics or groups 
of harmonics of these models, and of the uncertainty state 


of each parameter under consideration due to the model 
determination features. 

The comparison of model data with land point data on 
gravity anomalies and quasigeoid heights on the territory of 
Novosibirsk region is presented. 

On the territory under study 17 second-order gravity base 
network stations were chosen with gravity determination 
accuracy being + 0.05 mGal. The territory under investiga- 
tion includes 208 stations where geometric leveling (first to 
fourth order) was conducted and the values of normal heights 
were obtained. At the same points, satellite coordinates were 
determined for developing active base stations networks. As 
a result of satellite network adjustment geodetic heights were 
obtained, with mean square errors ranging from 1.5 cm to 
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Fig. 3 The results of calculating the resolution capabilities by the geopotential ultrahigh-degree models characteristics with approximation error 


of 1 cm and 1 mGal (period 2008-2019) 


3.1 cm, 1.8 cm, on average (Karpik et al. 2010). On the 
territory under consideration gravity anomalies variations 
and elevation differences are insignificant. 

The scheme of differences between the quasigeoid heights 
obtained from model EIGEN-6C4 and those of geoimet- 
ric leveling and GNSS measurements is given in Fig. 4. 
Applying model EIGEN-6C4 made it possible to achieve 
better results on the territory under consideration, with max- 
imum difference being 22 cm and mean square error about 
8 cm. 

The scheme of differences between the gravity anomalies 
obtained by EIGEN-6C4 model and those based on the 
ground data is shown in Fig. 5. 


Maximum difference amounted to 8mGal, with mean 
square error being about 4 mGal. 

The comparison of model and ground data for the gravity 
and quasigeoid heights anomalies on the territory of Novosi- 
birsk region (for 9 geopotential models) is given in Table 1. 


2 Conclusions. Interpretation 


Comparison of degree dispersions of geopotential harmonic 
coefficients and their errors is presented for 70 current 
global models built by the results of measurements of space 
gravimetric missions CHAMP, GRACE and GOCE. 
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Fig. 4 Scheme of differences between the heights of the quasigeoid obtained from the EIGEN-6C4 model and those obtained from geometric 


leveling and GNSS measurements, m 


The models under study revealed: the dependences of 
gravity anomalies and quasigeoid heights approximation 
errors on the degree of geopotential expansion into Fourier 
transform; ratio errors of geopotential harmonic coefficients 
determination; coefficients dispersions and those of gravity 
anomalies coefficient errors. 

Satellite models of the Earth’s gravity field built by only 
space gravimetry missions are not sufficient for obtaining 
quasigeoid heights with accuracy of 1-2 cm and spatial 
resolution less than 100 km. Satellite models are to be used 
in combination with land and topographic data to reestablish 
the global quasigeoid. 

Models resolution at the satellite altitude level satisfies the 
requirements on the models accuracy stated by the develop- 
ers. In case of necessity for recomputation of the gravity field 
characteristics for the Earth’s surface the specified resolution 
capability is reduced. 

The compared values of the model and ground data on 
the Novosibirsk region areas demonstrated that the least 
values of differences were obtained on the global geopo- 
tential models based on the GOCE mission data (Kanushin 
et al. 2015). The global models accuracy values for the 
territory of Novosibirsk region do not come up to 1 mGal 
and 1 cm, being in the range of 4-10 mGal and 7-20 cm, 
respectfully. 

The obtained accuracy values of the gravity anomalies 
and quasigeoid heights in current global geopotential models 


under study may be considered as maximum for the newly 
constructed models by the data of GOCE project. 

Space gravity mission GOCE made it possible to essen- 
tially improve harmonic coefficients for degrees 100-250. 

Harmonic coefficients of geopotential obtained for current 
EGF models are in conformity with each other within mean 
square errors. 

Further research of the global models is required using the 
results of space gravity mission GOCE as the most successful 
mission for investigating EGF. 

Accuracy characteristics of the models under study 
approximate the specified ones in process of their 
construction, but fall short of them as concerns the required 
resolution capability. 

Current global models of gravity field create new oppor- 
tunities and essentially extend the range of challenges in 
studying the Earth’s gravity field and its figure. 

Comparing the evolution of resolution capability and 
approximation accuracy as regards current global models of 
the Earth’s gravity field we may come to the conclusion that 
at the current stage of technological and methodical level, 
EGF model resolution capability and accuracy improvement 
has reached its maximum. For further improvement of EGF 
approximation accuracy new technological approaches and 
methodological principles should be developed. This implies 
that a new concept of the Earth’s gravity field study should 
be developed. 


Studying the Evolution of Resolution Capabilities and Approximation Accuracy of Global Models by Spectral. . . 113 


ve i 78° 79° 80° 81° 82° 83° 


Isa -3 


53 53° 


7 78° 79° 80° 81° 82° 83° mGal 


Fig. 5 Scheme of differences between the EIGEN-6C4 model reconstructed and ground values of gravity anomalies obtained for the Novosibirsk 
region in mGal 
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Abstract 


The paper reviews models of tidal and non-tidal variations of the Earth’s gravitational 
field. Proposing an algorithm for the estimation of the Stokes coefficients based on 
inter-satellite measurements of low-orbit spacecrafts. By processing measurements of the 
GRACE mission, we obtained experimental estimates of gravity field monthly variations. 
The analysis of these values was carried out by calculating the change in the equivalent 


water height for a given area. 


Keywords 


GRACE - Gravity field - Atmospheric tides 


1 Introduction 


In 2002 satellites of the GRACE mission (Greicius 2013) 
were launched into near-earth orbit, allowing the mapping 
of the Earth’s gravitational field at monthly intervals. The 
spatial resolution of derived models is several hundred kilo- 
meters, determining the height of a quasigeoid with an 
accuracy of a few centimeters. 

The basis for calculating the monthly models of the 
gravity field is the joint processing of inter-satellite mea- 
surements and data on the position of spacecraft obtained 
from measurements of the navigation receiver (Beutler et al. 
2010). The procedure consists of three stages: 

e for each spacecraft (SC) the trajectory of its motion 
obtained from measurements of the onboard receiver is 
aligned by a model; 

e asystem of equations is formed, linearized in the vicinity 
of the obtained orbits using inter-satellite velocities and 
orbits of both spacecrafts as measurements; 

e a joint system of linear equations is formed on the 
monthly interval and is solved by the least squares 
method. 
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As a result of the estimation, independent monthly models 
of geopotential are obtained. 

As is known, in addition to the tidal and relativistic 
effects corresponding to the IERS 2010 convention (Petit and 
Luzum 2010), the processing of spacecraft in low orbits takes 
into account the influence of non-tidal atmospheric loading 
calculated according to the NWP (Numerical Weather Pre- 
diction model) of the European Center for Medium-Term 
weather forecasts (ECMWF -— European Center for Medium- 
Range Weather Forecasts) at the German Research Center for 
Geoscience (GFZ German Research Center for Geoscience) 
(Dobslaw et al. 2016). 

In current work we compare the first results of in-house 
estimation of monthly gravity field models based on satellite- 
to-satellite measurements with other solutions, and also 
study the contribution of non-tidal atmospheric loading 
using a specially developed experimental programming and 
mathematical software KAGOR (KArtography Geodesy and 
Orbit determination). KAGOR allows simulating orbits, 
kinematic trajectories and corresponding inter-satellite 
measurements with any given parameters, as well as to 
process real satellite-to-satellite measurements to estimate 
gravity field models. 
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2 Non-tidal Atmospheric Loading 


Based on the analysis and forecasting of the data from 

the operational high-resolution global model of the numer- 

ical weather forecast of the European Center for Medium- 

Range Weather Forecasts and bottom pressure data obtained 

by simulation without imposing additional restrictions on 

the global general circulation model MPIOM and ECMWF 
data, a product of non-tidal gravitational variations in the 
atmosphere and oceans is formed (Atmosphere and Ocean 

De-Aliasing Level-1B, AOD1B) which provides a priori 

information on temporal variations in the gravity field due to 

global mass redistribution in the atmosphere and the ocean 

(Dobslaw et al. 2017). 

AOD 1B was determined by the fully normalized Stokes 
coefficients of the anomalous external gravitational field of 
the Earth caused by the redistribution of masses predicted by 
the numerical models described above. 

AOD 1B consists of four sets of coefficients: 

e ATM -— the effect of the atmosphere, which includes 
contribution of the atmospheric surface pressure over the 
continents, the static contribution of the atmospheric pres- 
sure above the surface of the oceans to bottom pressure 
and the much weaker contribution of density anomalies of 
the upper atmosphere over continents and oceans; 

e OCN -contribution of the movement of ocean masses to 
bottom pressure; 

¢ GLO -the sum of ATM and OCN coefficients; 

e OBA — model bottom pressure data that include contribu- 
tion of water and air over the entire surface of the Earth. 
This product is set with a resolution of 3 h. Since Jan- 

uary 1, 1976, atmospheric anomalies and bottom pressure 
anomalies have been arranged in a series of spherical func- 
tions up to the 100 d/o to ensure the processing of space 
geodetic missions equipped with laser retroreflectors, in 
particular, LAGEOS. Since January 1, 2000 AOD 1B is set 
by decomposition to a degree and order of 180 and, in 
addition to atmospheric anomalies and anomalies in bottom 
pressure, takes into account density anomalies of the upper 
atmosphere, which allows to meet the high requirements of 
modern space geodetic missions aimed at the study of gravity 
field (CHAMP, GRACE, GOCE). 

The influence of the atmosphere consists of two parts: 
non-tidal and tidal. Depending on the direct gravitational 
attraction of atmospheric masses by bodies generating tidal 
potential and, more significantly, on the periodic motion of 
the lower boundary of the atmosphere caused by the tide 
in the Earth’s solid body and oceans and on the absorption 
of solar radiation by the atmosphere, which causes a signal 
with a period close to 24 h. The tidal part includes 12 waves 
(P1, S1, K1, N2, M2, L2, T2, S2, R2, T3, S3, R3) and is 
calculated by updating the coefficients of the corresponding 
waves of the tidal potential in the observation interval 2007— 
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2014. A non-tidal atmospheric load is formed by subtracting 
the obtained waves from the measurements, the study of the 
influence of which this work is devoted to. 

An example of non-tidal atmospheric load for the Ist 
January, 2003, 00:00 in the form of an equivalent water 
height (EWH) Using a Gaussian filter (300 km) (Wahr 2007) 
is shown in Fig. 1. 


3 Estimation of Monthly Gravity Filed 
Models 


As previously noted, the gravity field monthly model estima- 

tion procedure consists of three stages: 

e for each spacecraft from a pair, its trajectory obtained 
from the measurements of the onboard receiver using 
the kinematic method of high-precision determination 
(Precise Point Positioning — PPP) is independently aligned 
by a model at a 3-h interval; 

e at the same time interval a system of equations is gener- 
ated, linearized in the vicinity of the obtained orbits using 
inter-satellite velocities and orbits of both spacecraft as 
measurements; 

e joint system of linear equations on a monthly interval 
being formed and solved by the method of least squares. 
Independent monthly models of geopotential are the result 

of an estimation. 

Any processing is based on a measurement model. In 
general, it has the following form: 

y(t) = f (t, ro, vo, Cnm, Sum, P,- -.) + €, d) 
where y(t) measurements, ¢ time, ro, vo initial conditions 
in inertial reference frame, Cam, Snm Stokes coefficients of 
degree and order less or equal then n and m, p instrumental 
parameters, € measurement error 

At the first stage, the kinematic trajectory of the spacecraft 
(the position of the spacecraft for all navigation receiver 
measurements epochs on the processing interval calculated 
with PPP technique) was used as measurements. In the sec- 
ond, besides the trajectories of both spacecraft, the relative 
velocity was added. 

At the first stage, we estimated the initial conditions of the 
spacecraft and the parameters of the accelerometer, while the 
Stokes coefficients was not estimated. For the calculation of 
the spacecraft coordinate system, star camera data was used. 
The orbit parametrization is generally similar to that used in 
the work devoted to the analysis of high-precision methods 
of ballistic-navigational support of space geodetic complexes 
using the Geo-IK2 spacecraft as an example (Zaliznyuk et 
al. 2019). At the second stage, Stokes coefficients up to the 
80th degree and order and model empirical parameters of 
the measurement errors of the inter-satellite line were added 
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Fig. 1 Non-tidal atmospheric load in the form of equivalent water level for the 1st January, 2003 00:00 


to the refined parameters. At both stages, the relationship 
between the dependent variables and the independent ones 
is nonlinear, so linearization is performed in the vicinity 
of the a priori orbit. The a priori orbit was calculated by 
numerical integration of the system of equations of motion 
with the initial conditions, taking into account a number of 
perturbations in accordance with the IERS 2010 convention. 
Table 1 shows the perturbations taken into account. At both 
stages, the derivatives of the measurements with respect to 
the orbital parameters were calculated by numerical integra- 
tion of the variations equations. 


Table 1 The perturbations taken into account when integrating the 
system of equations of motion 


Perturbation Model 
Static Earth gravity field EGM2008 
The direct attraction of the Sun, the DE430 


Moon and planets 


Solid Earth tides JERS Conventions 2010 
Ocean Loading FES2004 
Pole tide JERS Conventions 2010 


Deformation due to polar motion IERS Conventions 2010 
IERS Conventions 2010 
AOD1B RL06 


Accelerometer data 


Relativistic effects 
Non-tidal atmospheric loading 
Non-gravitational perturbations 


The contribution of various gravitational effects and a 
typical change in gravity filed over a monthly interval are 
shown in Fig. 2 in the values of the root of the degree variance 
in the logarithmic scales. Obviously, the static gravitational 
field has the greatest influence on the motion of spacecraft 
in the entire gravity spectrum. Then comes the ocean tide, 
whose contribution to the low-wave part of the spectrum is 
smaller than the changes in the gravity at monthly intervals. 
The atmospheric non-tidal load, for its part, has a greater 
effect on the gravity filed than its monthly variability up to 
10-15 harmonics and less on higher harmonics. 

Also, notice that when integrating the equations of 
motion, in addition to the conventional perturbations, 
measurements of non-conservative forces made by 
accelerometers installed on board both GRACE spacecraft 
were taken into account. Because of that, we need to add 
to the vector of the estimated parameters a systematic error 
for each of the three axes of the accelerometer and scaling 
factors in all of them as well. 

At the third stage, from the matrices of conditional equa- 
tions formed at the second stage corresponding to 3-h inter- 
vals, the matrix of normal equations was formed on a 
monthly interval. The vector of specified parameters of the 
monthly interval is a superposition of 3-h vectors in terms of 
orbital and instrumental parameters, and includes one set of 
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Fig. 2 The contribution of 
various gravitational effects 
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Contribution of various gravitational effects 
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Stokes coefficients for the entire interval. On average, one 
monthly solution includes about 3,600,000 measurements, 
10,000 estimated parameters, takes about 200 GB on the hard 
drive and lasts 3 days on PC. 


4 Results 


To assess the contribution of non-tidal atmospheric load to 
the estimation of monthly gravity field models, two series 
of solutions were carried out: with and without AODIB. 
In Fig.3 is a diagram of the solutions made to estimate 
the monthly gravity models. The missing of 2 months of 
both sets of solutions is due to the lack of accelerometer 
data at these intervals. This interval was chosen because of 
the relatively more successful functioning of the mission as 
a whole (a relatively small number of hardware failures). 
Considering the informativeness of the results of comparing 
the obtained time series of monthly gravity models and 


Fig. 3 Diagram of monthly 
gravity field solutions 


the high cost of computational resources for calculations, 
solutions without the use of non-tidal atmospheric loading 
for 2005 were not implemented. 

To evaluate the characteristics of both obtained time 
series of the gravity models, the equivalent water height 
was calculated using a Gaussian filter (300 km) in the 
Greenland region, since for this region a pronounced annual 
signal of changes in the gravity field and a long-period 
trend are observed. A rectangle was used (60..85 deg. of 
latitude; —70..—20 deg. of longitude). For the entire 3- 
year interval, the average value of gravitational potential 
was found, the pairwise difference of all monthly models 
with the obtained average was constructed, then for each 
point in the region with an increment of 1° the equivalent 
water level of the difference of the monthly model with 
the average was considered and the average value for the 
entire region was considered. The values obtained for both 
series of gravity models and similar values calculated from 
models of other official GRACE measurement processing 
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Fig. 4 The equivalent water level in Greenland (60..85 deg. of latitude; 
—70..—20 deg. of longitude), calculated and averaged over the grid 
with a step of 1 degree using a Gaussian filter (300 km) over a 3- 


centers (CSR, GFZ, JPL) are shown in Fig. 4. Considering 
that in the C2ọ series obtained from GRACE measurements, 
an unknown 161-day period is observed, the values of the 
C29 coefficient were taken from the results of processing 
five satellites (LAGEOS-1/2, Stella, Starlette, AJISATD) laser 
location data (Satellite Laser Ranging, SLR) (Loomis et al. 
2019). 

The results of the three official GRACE data centers are 
in good agreement with each other, due to the similarity of 
the models and processing algorithms used. The solutions we 
obtained taking into account the non-tidal atmospheric load- 
ing agree much better than the solutions obtained without 
taking into account AOD 1B. It can be seen that ignoring the 
non-tidal loading in processing introduces a distortion in the 
estimation of the gravity field models the same size as the 
signal of its change. This is the first result for the KAGOR 
software and this software has potential capability to increase 
the accuracy of the solutions obtained, in particular, in outlier 
search algorithms and reducing the influence of anomalous 
measurements. Nevertheless, this SW already provides an 
account of all the necessary effects and allows to obtain an 
adequate estimate of the gravity field signals. 


5 Conclusions 


1. The KAGOR (cartography, geodesic research and orbit 
determination) experimental software has been devel- 
oped, which allows to model orbits, kinematic trajectories 
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year interval according to monthly models of gravity filed by different 
analysis centers CSR, GFZ, JPL and IAC 


and to take corresponding satellite-to-satellite measure- 
ments with any given parameters, as well as to process 
real satellite-to-satellite measurements to estimate gravity 
field models. 

2. A general form of a measurement model was presented, 
the perturbations taken into account, the instrumental 
corrections taken into account, and the monthly solution 
generation algorithm. 

3. The contribution of non-tidal atmospheric loads to the 
monthly gravity filed models were analyzed. 

4. The monthly gravity field models were calculated taking 
into account and w/o taking into account non-tidal atmo- 
spheric loads. 

5. The results obtained show that the signal of non-tidal 
loads exceeds the signal of monthly variations of the 
gravity field. 
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Abstract 


The actual problem of creation of gravity maps which can be used as correctors in map- 
matching- navigation systems (MMNS) is considered. It is assumed to use information from 
two sources: available gravity maps and measurement data of various parameters of the 
Earth’s gravity field (EGF) (acceleration of gravity, deflection of the vertical, horizontal 
gradients of the potential). 

Obtaining the required amount of measurement data of a specified accuracy is considered 
as a problem of planning primary measurements. At the first stage of its solution, the 
number of measurements necessary to achieve the required accuracy of the parameters of 
the gravitational map is minimized. It takes into account the accuracy and performance of 
measuring devices. The second stage assumes that the human operator specifies the position 
of the points at which to measure the parameters of the EGF with a known accuracy. 

Joint processing of information is carried out on the basis of the collocation method. 
The estimation of EGF navigation quality is considered in a specified point and its nearest 
neighborhood. This approach is reasonable, if in the future it is required to optimally form 
the trajectory of the MMNS in a specified area or to a specified destination point. It is 
assumed rectilinear motion of the object from this point with a known initial level of 
coordinate errors. Position errors are defined as functions of distance from the initial point 
and direction of movement. 
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— instability of EGF parameters caused by various external 


1 Introduction factors, characterized by the value ~0.2 mGal (Nepok- 
lonov et al. 1996), and can be compensated; 

Currently, there is a growing interest in the use of EGF for -— recently, significant progress has been made in the devel- 

autonomous navigation of mobile objects. This is due to the opment and production of precise EGF component sen- 

following circumstances: sors (gravity, deflection of the vertical, horizontal gra- 

— gravitational field, as well as magnetic, has a global dients of the potential), functioning on mobile objects 
character, including underwater and outer space; (Nepoklonov et al. 1996; Peshekhonov et al. 2017). 

— EGF parameters, unlike magnetic ones, are practically not One of the actual problems of autonomous navigation with 
affected by interference (Beloglazov et al. 1985); the use of EGF is the task of gravity maps formation (Bel- 
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oglazov et al. 1985; Berkovich et al. 2018; Dzhandzhgava 
et al. 2013; Kotov et al. 2017). The basic requirements for 
such maps are determined by the level of errors, ensuring 
their effective use in the MMNS as correctors, as well as the 
initial alignment. Existing EGF maps and models do not fully 
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meet the accuracy requirements. Requirements are increased 
strongly due to the improvement of EGF parameter sensors. 

Precise measurement of EGF parameters (an order of 
magnitude more accurate than maps) at separately located 
points remains very costly at the present time. That is why it 
is important to improve methods of information processing, 
allowing to form (quickly update existing) gravity maps with 
the attraction of additional precise measurements. 

The article deals with three relatively independent prob- 
lems of ensuring the creation of gravitational maps for the 
required accuracy of MMNS. The first one is to ensure 
the specified accuracy of maps by performing minimum 
measurements of the gravitational parameters at sites on 
the Earth’s surface. The problem of optimal estimation of 
map errors using additional measurements is the basis of the 
proposed method for creating updated maps. Its peculiarity 
is the relatively small number of measurements performed 
mainly at sites on highways, as well as the different accuracy 
of measurements and the correlation of their errors. The 
problem of evaluating the navigation quality (informative- 
ness) of the navigation field is solved for the case when its 
results are used in forming the trajectories of the MMNS 
movement in a specified area or to a predetermined endpoint. 


2 Optimization of Projects for Primary 
Measurements in the Creation 
of Maps 


This optimization is due to the need to increase the accuracy 
of the parameters of the EGF with limitations on the possibil- 
ity of carrying out measurements at points called basic (BP) 
on the ground surface. These capabilities are established by 
the number of BP and the accuracy of the EGF parameters. 
In the BP the composition of primary measurements of 
EGF parameters is: altitude, acceleration of gravity, and 
components of the plumb deviation, horizontal potential 
gradients. As arule, in an arbitrary BP on the earth’s surface, 
it is assumed to measure the height and gravity. In addition, 
deflection of the vertical and horizontal gradients can be 
measured. The set of the parameters and the requirements 
to the accuracy of their measurement in a particular BP 
are the result of solving the optimization problem under 
consideration. BP’s position is a priori unknown. Therefore, 
an integrated index of planning quality is introduced, namely, 
the mean value of the root-mean square error (RMSE) of the 
EGF parameter in the local area of the terrain. It is used as an 
optimality criterion to be minimized. 

Among the limitations in the problem are the edges of 
the area, the number of the same type of BP (such as where 
the same set of the measured parameters of the EGF). The 
problem is to find the coordinates of all BP, for which the 
integrated quality index is minimal. The main difficulty of 
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finding its optimal solution is seen in a large number of 
optimized parameters — BP coordinates. It is also necessary 
to take into account that such solution is only the “first 
approximation” of the optimal solution of the problem. 
The final solution made by the human operator must take 
into account other significant factors that cannot or should 
not be included in the problem to be solved. Given this, 
it is permissible to search for a suboptimal solution. The 
proposed solution to the problem of BP placement, which 
provides a demanded accuracy of EGF parameters, includes 
two stages: 

— determination of the minimum number of BP, potentially 
providing a demanded accuracy of the parameters within 
a given area of the terrain; 

— specification of the position of each of the BP in a given 
area (“placement” BP). This stage is performed by the 
operator using cartographic materials (Sholokhov et al. 
2014). 

The minimum required number of BP is found through the 

density of the location of the BP on the ground surface, which 

is determined by the average distance between the BP. In 

General RMSE o, of geodesic parameter p or integral index 

can be represented by the function 

Op = fp (T, Om, Q), (1) 

where r — grid step, Om — RMSE of measurements in BP, Q — 
other cartographic source data BP coordinates, parameters 
of the covariance function of the gravity anomaly error. The 
function f, is formed using the RMSE o, of the EGF param- 
eters obtained by the collocation method for the values of its 
arguments specified in the grid nodes. Between nodes values 
Op are found by interpolation.The relationship between the 
desired RMSE o,, the number of BP, which is funded by the 
grid step r, and accuracy of additional measurements in BP 
allows to solve the inverse problem of finding the function’s 
argument, i.e. when o, is given. The program allows quickly 
solving the problem of planning the finding of these data, 
based on the specified levels of error (or the integral index of 
accuracy). The result of solving such problems provides the 
final specification of the BP position (their placement on the 
map) by the operator on the local terrain area. 


3 Methods for Solving the Problem 
of Creating Maps 


The creation of gravity navigation maps includes the solution 

of the following particular problems: 

— formation of sets of cartographic and newly measured 
gravimetric parameters; 

— estimation of map parameter errors in the areas where 
additional measurements were made; 
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Fig. 1 The example of gravimetric map (Bouguer anomaly) with additional measurement points 


— formation of an electronic map of the required gravimetric 
parameter; 
— evaluation of the navigation quality of the generated map 

(Kotov et al. 2017). 

Joint data processing of available EGF maps and precise 
additional measurements is based on the collocation method 
(Berkovich et al. 2018; Sholokhov and Druzhinin 2011; 
Moritz 1980). It allows finding the optimal estimates of the 
EGF parameters at an arbitrary point, taking into account 
additional measurements made beforehand at BP within the 
local area. Characteristics of the accuracy of the obtained 
estimates can also be found. 

One of the problematic issues of finding parameter esti- 
mates is the choice of a priori values of parameters charac- 
terizing the errors of EGF maps. It significantly affects the 
final estimates of the EGF parameters, since the amount of 
measurement data is small. To solve this problem, we use 
a special approach that allows us to take into account the 
optimal covariance of estimates of the desired navigation 
parameter at a given point, obtained with different combi- 
nations of values of a priori data on the errors of maps. 

The solution to this problem is realized by specialized 
software. A sample main program window depicting 
Bouguer anomaly layer with points in which additional 
measurements are made (marked with the symbol A), shown 
in Fig. 1. 

Program in conjunction with geographic information sys- 
tems (Professional GIS “Panorama” 2019) provides a solu- 
tion to the following applied tasks: import and visualization 
of gravimetric information, manual input of information 
from the various sensors, the joint processing of the mea- 


surement results, the formation of databases with the visual- 
ization of the results, estimation of accuracy of geodetic and 
gravimetric parameters, formation navigation maps layers of 
different gravimetric parameters (Bouguer and free-air grav- 
ity anomaly, the components of the deviation of the plumb 
line, the second derivative of the potential), adjustment of 
the maps according new measurement data, the reduction 
of gravimetric parameters at points with different heights to 
determine a gravimetric parameter of an arbitrary point of the 
route at the specified values of coordinates, evaluation of the 
navigation quality of the navigation field. 


4 Estimation of the Navigation Quality 
of the Navigation Field 


The informativeness of the EGF parameter is determined at 
a specified point, from which the rectilinear motion of the 
object in a known direction is assumed. This approach differs 
from the traditional one (Beloglazov et al. 1985; Koneshov 
et al. 2016; Peshekhonov et al. 2017). It allows you to 
prepare the initial data for the subsequent optimal formation 
of the trajectory of the MMNS in a specified area or to a 
specified endpoint. An indicator of navigation quality is the 
RMSE of the object’s coordinates. Finding simple analytical 
expressions that allow to obtain achievable levels of RMSE 
of coordinates depending on these factors is difficult and in 
practice impractical. 

Therefore, a numerical assessment of these levels is car- 
ried out, based on modeling the process of correction the 
coordinates of the MMNS. In it random values of the data 
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Fig. 2 Example of navigation quality estimation for the specified point 


are generated in accordance with the required RMSE. You 
need to find the following functions fx, fy for RMSE ox, oy 
of coordinates. They are allow us to solve inverse problems — 
to find arguments of functions when levels of RMSE of coor- 
dinates are given. Such problems are of the principal interest, 
since in practice the levels of RMSE of required parameters — 
coordinates are set, and the conditions for achieving the 
desired accuracy, determined by the arguments of functions, 
are subject to definition. The solution of inverse problems is 
obtained by finding the roots of nonlinear equations 


ox — fx (Og, Tg, Om, Le)| a= 
oy — fy (Og, Tg, Om, Le) ogr = 


l (2) 


where Og, rg — RMSE and spatial correlation radius of 
the EGF errors are specified, the accuracy Om of the EGF 
parameters measurements is known, and the length L, of the 
MMNS correction section is to be determined. The result of 
navigation quality estimation is presented in the form of a 
map of the minimum possible (achievable) location error for 
one point of a given object trajectory. An example is shown 
in Fig. 2. The contour lines show the achievable levels of 


distance RMSE o = ,/o% + 07 in the rectilinear motion of 
the object from the starting point of the trajectory (marked 
“+”) in the appropriate direction. The error levels of the 
coordinates at the reference point, as well as other parameters 
of the object motion and the parameters of the MMNS are set 
by the operator. 
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5 Conclusion 


The problems of creating gravity maps to specify the position 
of movable objects are considered. It is planned to use 
information from two sources: available gravity maps and 
measurement data of various parameters of the EGF. Joint 
data processing is performed on the basis of the standard 
collocation method. The peculiarity of the application of this 
method to find estimates of the EGF parameters is revealed. 
The estimation of EGF navigation quality is considered in 
relation to a given point and its nearest neighborhood. This 
approach is reasonable, if in the future it is required to 
optimally form the trajectory of the movable object in a 
specified area or to a given endpoint. 
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Abstract 


The current state of satellite altimetry development allows the use of altimetry data in 
determining the detailed characteristics of the Earth’s gravity field on the surface of the 
ocean in the form of models of geoid heights, deflection of vertical and gravity anomalies. 
The article presents the bistatic altimetry method based on the reflected signals of global 
navigation satellite systems (GNSS). In this method, measurement redundancy is necessary 
to solving the problem of determining the height of the geoid with high spatial resolution. 
The experiments showed the possibility of using the code and phase of the received signal. 
The experiments showed the possibility of using method of bistatic GNSS-altimetry to 


determine the height of the geoid. 


Keywords 


Altimetry - Geoid - GNSS-R - Measurement methods 


1 Introduction 


The definition of the geoid is one of the main tasks of 
geodesy, gravimetry and oceanography. Geoid is a surface of 
equal gravitational potential on the Earth, containing a point 
taken as the beginning of the height count. The surface of the 
geoid coincides with the surface of the World Oceans in the 
absence of disturbing forces such as wind, ocean currents, 
tides and conditionally continues under the continents. The 
height of the geoid is the elevation of the geoid above the 
ellipsoid. 

Currently, the method of active satellite altimetry is used 
to measure the height of the geoid on the oceans. This method 
is based on the use of radio altimeters placed on board special 
spacecraft for geodetic purposes. 

Satellite altimetry measurements with improved accuracy 
and spatial resolution can come closer in accuracy and 
resolution to the results of marine gravity surveys. 
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2 Ways to Improve the Accuracy 
of Altimetry 


To date, over 10 national and international satellite altimetry 
projects have been implemented. Foreign altimetry missions 
are shown in Table 1. 

To improve the spatial resolution of altimetry measure- 
ments in the medium term, it is necessary to use: 


e altimeter with one carrier frequency in the Ka-band; 
e aperture radar synthesis method; 

¢ interferometric method; 

e constellation of spacecraft; 

e method of bistatic GNSS-altimetry. 


Satellite altimetry based on active monolocators has high 
measurement accuracy. However, the following disadvan- 
tages of active satellite altimetry can be identified: 


e the onboard radio altimeter consumes more power; 
e only one observation track of the current altitude profile 
in one pass of the spacecraft is realized. 


It is suggested that GNSS-R be used to eliminate the dis- 
advantages of active altimetry. GNSS-R is a remote sensing 
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Table 1 Existing and future foreign altimetry missions 
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Mission Height, km Inclination Accuracy, m Review period, days 
GEOS-3 (04.75-12.78) 845/115 0.5 - 
Seasat (06.78—09.78) 800km 800/108 0.1 3 
Geosat (03.85-09.89) 785.5/108.1 0.1 17.5 
ERS-1(07.91—03.2000)ERS-2(04.95—09.2011) 785km 785/98.5 0.05 35 
TOPEX/Poseidon(09.92—10.05)Jason-1 (12.01—06.13) 1336/66 0.025 9.92 
GFO (02.98-09.08) 800/108 0.1 17.5 
Envisat (03.02—04. 12) 785/98.5 0.03 35 
Jason-2 (06.08 — present) 1336/66 0.025 9.92 
CryoSat-2 (04.10 — present) 717/92 0.05 369 
HY-2A (08.11 — present)HY-2B (10.18 — present) 971/99.3 0.02 14,168 
SARALA/AItiKa (02.13 — present) 785/98.5 0.02 35 
Jason-3 (01.16 — present) 1336/66 0.01 9.92 
Sentinel-3A (02.16 — present)Sentinel-3B (04.18 — present) 814.5/98.65 0.02 27 
CFOSAT (10.2018 — present) 500/90 - 13 


method of the Earth that uses GNSS signals reflected from 
the water surface, allowing simultaneous observations at 
several points in a wide band. This method is promising for 
mesoscale altimetry. GNSS altimetry advantages: 


e many simultaneously processed altitudes (number of visi- 
ble GNSS satellites 20-30) (Sahno et al. 2009); 

e low power consumption and low weight of on-board 
equipment; 

* opportunity to determination of the speed and direction of 
the near-surface wind (Awange 2018); 

e presence of algorithms for signal processing (Jin and 
Komjathy 2010). 


3 Operating Principle 


The GNSS altimetry bistatic system is based on the reception 
of direct signals from navigation satellites by an antenna 
directed at the zenith. Signals reflected from the ocean 
surface in the specular zone are received by an antenna 
directed to the “nadir”, The operating principle is presented 
in Fig. 1 (Zavorotny and Voronovich 2000; Gleason et al. 
2005). To solve the problem of determining the height to the 
reflecting surface, it is necessary to determine the difference 
in the propagation time of the direct and reflected signals, 
the coordinates and speeds of the receiver and the satellite. 
The accuracy of determining the height depends on the 
noise of the equipment, the accuracy of accounting for the 
ionospheric delay, tropospheric delay, processing algorithms, 
etc. (Camps et al. 2017; Li et al. 2016). 


4 Methods of Bistatic GNSS-Altimetry 


Today we can talk about four methods: 


(1) Code method. Code method can be called traditional. 
Receivers of bistatic altimetry system determine 


pseudorange and pseudodoppler frequency bias for 
direct and reflected signals, as well as commercial 
receivers. The received direct and reflected signals 
are correlated with a replica of the open code 
signal generated in the receiver. The method ignores 
the encrypted signal components, thereby reducing 
the achievable accuracy of bistatic GNSS altimetry 
(Zavorotny et al. 2014). 

(2) Interferometric method. An interference method can be 
used to improve measurement accuracy (Martin-Neira 
1993). It is based on mutual correlation processing 
of direct and reflected signals from one satellite, 
received by two antennas with a large gain. Such 
antennas are used for spatial selection of received 
signals; 

(3) Method of SNR. The SNR method analyzes the change 
in the ratio of the signal (sum of direct and reflected 
signals) power level to the noise power depending on 
the satellite elevation angle, the signal wavelength, and 
the height of the receiving antenna (Lofgren and Haas 
2014). Measurements by this method can be imple- 
mented using a single antenna to receive direct and 
reflected GNSS signals. The dependence of the SNR 
value on the height of the receiving antenna in accor- 
dance with (1) 


SNR(t) = Acos (Zn cos (0 (t)) + v) : (1) 


where 6(t)—elevation angle of the satellite at time t, A— 

signal wavelength, h—height to the reflecting surface, ~ o— 

initial phase, A—amplitude of the total signal. Unknown 

parameters are h, @ o, A. 

(4) Phase method. Phase measurements have higher accu- 
racy than code measurements, but such measurements 
are ambiguous due to the uncertainty of the initial integer 
number of phase shift periods (Semmling et al. 2014). 
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GLONASS/GPS/ 
Galileo 


Fig. 1 Principle of operation of bistatic GNSS-R altimeter 


Fig. 2 Place of experiment 
p.1—place of measurements with 
the code method; p.2—place of 
measurements with the SNR 
method 


However, if it is necessary to determine the change 
(increment) in height to the reflecting surface, the phase 
measurements can be used as more accurate on condition 
that there are no cycle slips. The change in height from 
time ¢ ; to t ; is proportional to the change in the 
pseudophase difference: 


__ Ag) __ A” (ti) (2) 
2sin (0 (t;)) 2sin(@(t;))’ 

where Ag—difference between the pseudophases of the 
direct and reflected signals.It is necessary that the measure- 
ments satisfy the Rayleigh criterion reference required. The 
Rayleigh criterion is fulfilled at low elevation angles of the 
spacecraft or at a small surface roughness. The accuracy of 
phase measurements can be centimeters. 


5 Experimental Verification 
of the Methods 


Experimental measurements were performed to test bistatic 
GNSS altimetry methods from a bridge over a reservoir. 
In the first case, the GNSS altimetry code method was 
tested. The place of measurements is shown in Fig. 2. The 
receiving antennas were located at point 1. The results of 
measurements of the height to the reflecting surface are 
preceded in Fig. 3. From the measurement results, we can 
conclude that the optimal elevation angles of satellites for 
this method are 30°—90°, and the measurement error may be 
less than 1 m. 

The SNR method was tested in measurements at p. 2, the 
real and measured height of which to the reflecting surface 
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was 17.1 m, Fig. 4. Such a noticeable difference between the 
heights of p.1 and p.2 is associated with the difference in 
the height of the bridge and different levels of the reservoir’ s 
water in different seasons in which measurements were 
taken. 

An anechoic shielded chamber was chosen for testing 
the phase method. The measurement scheme is shown 
in Fig. 5. Providing the necessary roughness of the 
reflecting surface and the necessary elevation angles 
of the satellite is not always a feasible task for field 
measurements. 

Semi-natural modeling using a GNSS simulator in ane- 
choic chamber conditions allows (Frolov et al. 2017): 


e to minimize reflections of radio signals from obscuring 
objects; 

e to form a spatial navigation field of GNSS signals; 

e to change the signal power level; 


30 
Heigth, m 


e perform measurements both on the code and phase delay 
of the signal. 


Antennas receiving direct and reflected signals were located 
on a mast with a variable height. GNSS satellites sig- 
nals were simulated by a simulator manufactured by Navis. 
The results of measuring the difference between the pseu- 
dophases of the direct and reflected signals with a decrease 
and a further increase in mast height are shown in Fig. 6. 
The standard deviation of measurements was no more than 
0.2 cm. 

The phase altimetry measurements with the CYGNSS 
bistatic GNSS system over two independent passages above 
Qinghai Lake are shown in (Li et al. 2018). As the water 
surface of the lake approaches the equipotential surface of 
gravity, altimetry measurements above the lake can give an 
independent estimate of the height of the geoid along the 
track. The measurement data are relative due to an unknown 
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Transmit antenna 


Receiving 
antennas. 


Fig. 5 The scheme of the experiment in an anechoic chamber 


Pseudophase difference 


Pseudophase difference, m 


Fig. 6 The scheme of the experiment in an anechoic chamber 


integer number of phases. The measurements obtained after 
introducing all the necessary corrections were compared 
with the EGM2008 model along the track. The measurement 
results well repeat EGM2008, but anomalies of + 10 cm are 
distinguished, which are repeated in time for different GNSS 
satellites. 


6 Summary and Conclusions 


Thus, the experiments show the consistency known methods 
of the bistatic GNSS-altimetry. An anechoic chamber and a 
GNSS signal simulator can be used to simulate the operation 
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of a GNSS altimetry system using any of four known meth- 
ods. The highest requirements for a GNSS signal simulator 
are imposed during phase measurements. The accuracy of 
determining the height can reach less than 0.1 m with space 
placement of a bistatic GNSS altimeter. 
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Abstract 


Gravity changes associated with volcanic processes occur over a wide range of time 
scales, from minutes to years and with magnitudes between a few and a few hundred 
microGal. High-precision instruments are needed to detect such small signals and both time- 
lapse surveys along networks of stations, and continuous measurements at single points, 
are accomplished. Continuous volcano gravimetry is mostly carried out through relative 
gravimeters, either superconducting instruments, providing higher quality data, or the more 
widely used spring meters. On the other hand, time-lapse surveys can be carried out with 
relative (spring) gravimeters, that measure gravity differences between pairs of stations, 
or by absolute gravimeters, capable of measuring the absolute value of the gravitational 
acceleration at the observation point. Here we present the state-of-the-art of terrestrial 
gravity measurements to monitor and study active volcanoes and the possibilities of new 
gravimeters that are under development. In particular, we present data from a mini array of 
three iGrav superconducting gravimeters (SGs) at Mount Etna (the first network of SGs ever 
installed on an active volcano). A comparison between continuous gravity measurements 
recorded through the iGrav#016 superconducting gravimeter at Serra La Nave station 
(1730 m a.s.l.) and absolute gravity data collected with the Microg LaCoste FG5#238 
gravimeter in the framework of repeated campaigns is also presented. Furthermore, we 
introduce the Horizon 2020 NEWTON-g project (New Tools for Terrain Gravimetry), 
funded under the FET-OPEN Research and Innovation Actions call, Work Programme 
2016-2017 (Grant Agreement No 801221). In the framework of this project, we aim to 
develop a field-compatible gravity imager, including an array of low-costs Micro-Electro- 
Mechanical Systems (MEMS)-based relative gravimeters, anchored on an absolute quantum 
gravimeter. After the design and production phases, the gravity imager will be field-tested 
at Mt. Etna (Italy) during the last 2 years of the project. 
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High-precision gravity measurements provide a powerful 
tool for volcano monitoring, since they highlight processes 
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In spite of its potential, volcano gravimetry is not widely 
adopted by research institutions and agencies in charge of 
monitoring active volcanoes. This depends on several factors, 
including: (1) the cost of instrumentation, (2) problems that 
are inherent with the use and deployment of instruments 
intended for laboratory conditions in the harsh environments 
that characterize the summit zones of most active volca- 
noes, and (3) difficulty in interpreting gravity changes that 
may result from volcanic processes, as well as hydrological 
effects (changes in the undergroung water mass), and instru- 
mental artefacts. These difficulties are not insurmountable 
and efforts in a variety of volcanic settings highlight the 
value of time-variable gravimetry for understanding volcanic 
hazards, as well as revealing fundamental insight into the 
volcano dynamics (Carbone et al. 2017). 

The Istituito Nazionale di Geofisica e Vulcanologia 
(INGV) has operated a relative gravity network for the 
monitoring of Etna volcano since 1986 (Budetta et al. 1989). 
The network has been developed and has evolved over 
the years and currently it consists of: (a) 71 benchmarks, 
covering an area of about 400 km’, for relative gravity 
campaigns (LaCoste & Romberg model D and Scintrex 
CG-3M/CG-5 gravimeters were used over time); (b) three 
continuously running gravity stations equipped with iGravs 
superconductive gravimeters by GWR Instruments, Inc.; 
(c) 14 stations for absolute gravity (AG) measurements 
using the Microg LaCoste FG5#238. In all stations for 
AG measurements and in some stations for relative 
measurements, the vertical gravity gradient is also measured. 

The whole Etna gravity network is routinely occupied 
every summer (in the same season of the year to minimize 
seasonal variation effects). More frequent measurements 
are carried out through relative spring gravimeters in same 
elements of the network (almost monthly along the East- 
West profile; Fig. 1; Carbone et al. 2009; Greco et al. 2010; 
Del Negro et al. 2013; Bonforte et al. 2017) or during 
volcano activities. Since 2007, several absolute field surveys 
have been routinely carried out at Etna volcano to check 
the long term gravity variations and to confirm the gravity 
changes obtained by relative measurements (Pistorio et al. 
2011; Greco et al. 2012, 2015). To collect the absolute 
gravity data, the IMGC-02 absolute gravimeter (D’ Agostino 
et al. 2008), developed by Istituto Nazionale di Ricerca 
Metrologica (INRiM, Italy), was used in the period 2007- 
2009. Since 2009 up to date the Micro-g FG5#238 absolute 
gravimeter is used. 

In this paper, we present data collected with SGs and 
absolute gravimeters at Mt. Etna. Furthermore, we introduce 
the field-compatible gravity imager, that will be developed 
under the NEWTON-g project, including an array of low- 
costs MEMS-based relative gravimeters, anchored on an 
absolute quantum gravimeter. This measuring system will 
provide continuous imaging of gravity changes associated 
with variations in subsurface fluid properties, with unpar- 
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alleled spatio-temporal resolution and, consequently, with 
fundamental implications for risk management. We will 
deploy the new gravity imager at Etna volcano (Italy), 
where frequent gravity fluctuations, easy access to the active 
structures and the presence of a multiparameter monitoring 
system ensure an excellent natural laboratory for testing the 
new tools. 


2 The Mini-Array of Three iGrav SGs 
at Mount Etna 


Continuous gravity measurements have been successfully 
carried out at a number of volcanoes around the world using 
spring gravimeters. Nevertheless, these instruments do not 
provide reliable continuous data over intervals of weeks or 
more, because they are influenced by environmental factors 
and are subject to instrumental drift. Accordingly, most stud- 
ies of continuous gravity at active volcanoes have focused on 
the analysis of changes over time-scales of minutes to a few 
days (Branca et al. 2003; Carbone et al. 2006, 2013, 2019). 

An alternative to spring gravimeters for continuous mea- 
surements is given by superconducting gravimeters (SGs) 
that feature a much higher precision and stability than spring 
gravimeters. SGs are free from instrumental effects and thus 
allow to track even small gravity changes (1-2 u Gal) over 
a wide range of time scales (minutes to months). However, 
even the most portable SGs (e.g., the iGrav’ by GWR) are 
not ideal for installation in the vicinity of active volcanic 
structures. Indeed, they require AC power at the installation 
site and some kind of hut or vault to house the instrumenta- 
tion (Carbone et al. 2019). 

At Mt. Etna, the installation of a mini-array of three SGs 
(distances of 3.5-15.5 km from the active craters; Fig. 1) 
began in September 2014, when the iGrav#16 was installed 
in the facilities of the Serra La Nave (SLN) Astrophysical 
Observatory (1,730 m above sea level (asl); 6.5 km from the 
summit craters; Fig. 1). In the summer of 2016, iGrav#20 and 
iGrav#25 were also installed at La Montagnola hut (MNT; 
2,600 m asl; 3.5 km SE of the summit craters; Fig. 1) and in 
the facilities of INGV — OE located in the village of Nicolosi 
(NIC; 720 m asl; 15 km SE of the summit craters; Fig. 1), 
respectively. 

To our knowledge, these are the first SGs ever installed on 
an active volcano. 

Signals from these instruments have shown hydrologi- 
cally-induced components superimposed on gravity changes 
that are related to volcanic processes (Carbone et al. 2019). 

Figure 2 shows the signals acquired in the period 23 July— 
05 September 2016 at NIC (top; iGrav#25), SLN (middle; 
iGrav#16) and MNT (bottom; iGrav#20) stations, after the 
corrections for the effect of Earth tides (Wenzel 1996), local 
atmospheric pressure variations (Merriam 1992), and polar 
motion (Hinderer et al. 2015), are accomplished. 
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Fig. 1 Sketch map of Mt. Etna showing the mini-array of three SGs 
has been in continuous operation at Mount Etna since the summer of 
2016. The stations belonging the East-West profile (blue circles) where 


The reduced signals show common anomalies over inter- 
vals of 10-15 days. In particular, a gravity decrease is 
observed in the signals from SLN and MNT stations during 
August 8-18. The MNT/SLN amplitude ratio during this 
period points to the activation of a mass source located a few 
kilometers below sea level. 

Figure 2 also shows positive variations on time scales of 
the order of a few hours on the signal acquired at MNT 
(red arrows). Since these variations (average amplitude equal 
to about 2 microGal) are not present on the SLN and 
NIC signals and occur during heavy rainfall, they are most 
probably due to local hydrological effects. 


3 Comparison Between the Time Series 
from a SG and Absolute Gravity Data 


Figure 3 shows a 56-month long time series from iGrav#16 
SG (at SLN station), after removing the effect of Earth 
Tides, atmospheric pressure and polar motion. Figure 3 also 
shows absolute gravity observations (red points) at the same 
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x 


more frequent measurements are carried out through relative spring 
gravimeters are also indicated 


station (the total error on AG measurements at this station is 
typically + 3 microGal). 

Once a linear fit, obtained through AG data, is removed, 
there is a good fit between the two data sets. Continuous 
gravity changes are within 1—2 microGal of gravity changes 
measured using AG at the same site. This comparison allows 
to estimate the long-term drift of the iGrav SG, which is of 
the order of 9 microGal/year (Fig. 3). 

This example shows how repeated AG measurements 
and continuous gravity observations can be used together 
to get a fuller and more accurate picture of the temporal 
characteristics of the studied processes. 


4 New Tools for Terrain Gravimetry 


Gravimetry is a powerful geophysical tool that, through 
sensing changes in subsurface mass, can supply unique 
information on the dynamics of underground fluids, like 
water, magma, hydrocarbons, etc. This is critically important 
for both resource management and risk reduction. In spite 
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Fig. 2 Gravity time series from NIC (top), SLN (middle) and MNT stations during the 23 July to 03 September 2016 interval, corrected for the 
effect of Earth tides and local atmospheric pressure. The time series are also high - pass filtered (cutoff of 0.01 Hz) 


of this potential, gravimetry is currently underexploited. 
Indeed, continuous gravity measurements have been rarely 
performed, owing to the fact that the available instruments 
are not well suited for continuous measurements under 
severe field conditions. In addition, the high cost of available 
gravimeters limits deployment to, at most, a few sensors in a 
given area, implying that the achievable spatial resolution is 
always lower than needed. 

To overcomethe current limits of terrain gravimetry, 
the NEWTON-g project, that received funding from the 
European Commission’s Horizon 2020 programme, under 
the FET-OPEN 2016-2017 call, proposes the development 
of a new generation of gravimeters, based on MEMS 
and quantum technologies (www.newton-g.eu). MEMS 
technology will make possible the production of a device 
that is two orders of magnitude cheaper and lighter 
than currently available gravimeters. In the framework 
of NEWTON-g, one goal is focused on strategies to 
increase the stability of ordinary MEMS accelerometers, 
thus making them capable of measuring, beside inertial, also 
gravitational acceleration, with a sensitivity of a few tens of 
microGal/hz!? (Middlemiss et al. 2016). 

The quantum gravimeter that will be developed under 
NEWTON-g will be the ideal complement to the MEMS 
devices. We expect that the new quantum gravimeter would: 
(a) measure the absolute value of the gravity acceleration 


(absolute instrument), with a sensitivity and a stability at the 
uGal level; (b) perform continuous measurements at a high 
rate (2 Hz) over a long time interval (several years). 

A network of low-cost MEMS gravimeters, anchored 
to the absolute quantum device, will form the so-called 
gravity imager, which will provide imaging of sub-surface 
mass changes with unparalleled spatio-temporal resolution 
(Fig. 4). 

Once developed, the gravity imager will be field-tested at 
Etna volcano (Italy) during the last 2 years of the NEWTON- 
g project (summer 2020 — summer 2022). Insights from the 
gravity imager will also be used for volcanic hazards anal- 
ysis, to demonstrate the importance of using gravity to face 
problems of societal relevance. A successful implementation 
of NEWTON-g will open new doors and represents a fun- 
damental step that can move gravimetry into a cornerstone 
resource for geophysical monitoring and research. 


5 Conclusions 


Many geophysical phenomena are associated with mass 
transport and may induce gravity changes measurable at the 
surface over a wide range of time scales. Gravimetry can thus 
supply unique information on mass transport within the Earth 
system. 
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Fig. 3 (Top) time series from iGrav#16 at SLN station during the 
September 2014—April 2019 interval. Data are corrected for the effect of 
Earth tides, local atmospheric pressure and polar motion (grey signal). 
The time series is also high — pass filtered (cutoff of 0.01 Hz). The green 


Through the experiments accomplished at Etna during 
the past few decades many steps forward have been made 
towards the regular acquisition of high-quality gravity data. 
The main issues with continuous gravity measurements at 
active volcanoes is to obtain and maintain data to a suit- 
able standard (as for quality and continuity) against an 
adverse environment (high altitude, inaccessibility for sev- 
eral months, lack of mains electricity for power, variable 
temperature, pressure and humidity, seismicity, corrosive 
gases) and using instruments which are intended for use 
under laboratory conditions. 

The combined use of discrete (absolute and relative) and 
continuous gravity measurements is a unique tool both for 
studying the internal dynamic of a volcano and for surveil- 
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curve is the same signal, after correction for a linear trend deduced from 
the AG data (red points). (Bottom) same signals as in the top panel, but 
during the January 2017—April 2019 interval 


lance purposes. Such combined system has already allowed 
important conclusions to be drawn at Mt. Etna. 

The state-of-the-art instruments are, in most cases, heavy, 
expensive, power-hungry and difficult to operate. This pre- 
vents the installation of dense networks of continuously 
running gravimeters, implying that suitable spatio-temporal 
resolution cannot be attained. 

Less expensive and easier to deploy instruments will 
open new horizons for terrain gravimetry. In particular, the 
gravity imager that will be developed and field-tested under 
the NEWTON-g project will demonstrate the possibilities 
of a new generation of MEMS and quantum gravimeters, 
allowing a step that can turn gravimetry into a cornerstone 
resource for volcano monitoring and hazard assessment. 
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MEMS gravimeter 


Possibility to deploy: (for the first 
time!) a dense array of continuously 
running gravimeters. 

High spatio-temporal resolution 
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Absolute quantum gravimeter 


Possibility to have an absolute 
reference that can be used for 
continuous and discrete measures. 

Moveable absolute reference 


Fig. 4 The NEWTON-g gravity imager (currently under development) consist of an array of low-cost MEMS relative gravimeters, anchored to 


an absolute quantum gravimeter 
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Abstract 


The paper analyzes the effect of preliminary processing of gravity measurements on 
the accuracy of the marine gravity-aided navigation. The preliminary processing of the 
measurements is implemented in the filtering and smoothing modes. Obtained results are 
illustrated by a one-dimensional example of gravity-aided navigation problem. 
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1 Introduction 


Map-aided navigation systems using geophysical fields for 
a position update are employed for a wide class of vehi- 
cles (Stepanov 1998; Bergman 1999). The vehicle position 
is updating through the comparison of the measured and 
reference samples (profiles) of the field along the vehicle 
trajectory. In marine navigation it is common to use the 
field of Earth gravity anomalies (GA), since it is stable in 
time and can be measured using well-developed inertial sen- 
sors: high-precision accelerometers and gravimeters (Bishop 
2002; Sokolov et al. 2019). 

However, the measurement of the gravitational field 
onboard the moving vehicles has its own peculiarities: 
the signal measured by the gravimeter consists of both 
GA and inertial accelerations of the object, considered as 
errors. There are two possible ways for solve the gravity- 
aiding problem: use raw measurements and consider the 
detailed error model in the map-aiding algorithm or use the 
preprocessed measured samples in which the GA along the 
vehicle trajectory is estimated some way. In the second case, 
often used in practice, the measurements of the gravimeter 
are pre-processed (Wu et al. 2017). After that, the amplitude 
of the measurement errors of the GA significantly decreases, 
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Preliminary 


which allows us to simplify the map-aiding algorithm. At the 
same time, as it is shown in (Nosov and Stepanov 2018), the 
preliminary processing of measurements and the subsequent 
simplification of the algorithm can lead to an increase in the 
navigation error. 

This paper analyzes the effect of preliminary processing 
of gravity measurements on the accuracy of the marine 
gravity-aided navigation. It continues the authors’ study on 
a similar problem for underwater terrain-aided navigation in 
a presence of a white noise measurement error (Nosov and 
Stepanov 2018). 


2 Optimal and Suboptimal Solutions 
of the Gravity-Aided Navigation 
Problem 


Following (Nosov and Stepanov 2018), we consider 
a gravity-aided navigation problem in the simplest 
formulation, namely, we need to estimate a constant scalar 
random variable A (for example, a navigation system 
(NS) error from one of the coordinates) using a previously 
constructed GA map. This problem can be stated in Bayesian 
framework as follows (Stepanov 1998; Bergman 1999): 
to estimate the unknown parameter Eq. (1) using p scalar 
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measurements Eq. (2) 


A; = Ai- = A, (1) 


vy =o(% —A)+e = (A+ 4, (2) 
where x; are coordinates, provided by the reference navi- 
gation system, e.g. inertial one, in a discrete set of points; 
£; are random measurement errors; i = Tack d(*) is a 
exactly known function (map) describing the dependence of 
the field on coordinate. For simplicity, we consider A and ¢; 
as Gaussian with known stochastic properties. 

To solve this problem in optimal way, the well-known 
numerical point-mass and sequential Monte-Carlo methods 
are used (Stepanov 1998; Bergman 1999; Nordlund 2002; 
Gustafsson et al. 2002). However complex behaviour of 
measurement errors £; which in turn requires high-dimension 
model makes the estimation algorithm implementation com- 
putationally expensive, because nonlinear Bayesian estima- 
tion methods are subject to “curse of dimensionality” (Daum 
and Huang 2003). In this case it is suitable to pre-process 
measurements to estimate the GA along the vehicle tra- 
jectory, i.e. the value of ¢;(A). The measurements Eq. (2) 
are represented as a sum of some random process which 
describes the field values g; = (x; — A) and the error 
same as above y; = g; + £i. The problem of estimating g; 
is linear one, which allows the use of Kalman-type filtering 
and smoothing algorithms (Kalman 1960; Gelb et al. 1974; 
Peshekhonov and Stepanov 2017). 

After the preliminary processing, the measurements used 
to solve the navigation problem can be written as follows: 

Îk = h(x — A) + Sk, (3) 
where jx is the estimate, provided by the preliminary filter- 
ing or smoothing algorithm; sx is the corresponding estima- 
tion error; k =1...n. 

Note that preliminary processing in itself does not make 
the algorithm for estimating A simpler. Moreover, if we use 
the whole set of measurements Y, after the pre-processing 
and take a proper account of the estimation errors, we can 
show that the estimation accuracy of A will remain at the 
same level. However, since preliminary processing signifi- 
cantly reduces the measurement errors it becomes possible 
to reducethe number of measurements used to estimate A 
from p to n. This is achieved by increasing the interval Atı 
for measurements Eq. (2) compared with the interval At for 
measurements Eq. (3). Furthermore if the interval between 
measurements Eq. (3) is chosen to exceed the correlation 
interval for the error sg, then its model can be approximated 
by discrete white noise, the level of which will correspond to 
the solution of a filtering problem or a smoothing one. This 
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simplifies the model used in nonlinear algorithm, since there 
is no need to describe the complex behavior of measurement 
errors. 

As indicated in the introduction, such a two-stage scheme 
for gravity-aided navigation algorithm can lead to an increase 
in positioning errors. To evaluate losses in accuracy, we will 
use the procedure based on comparing the unconditional 
calculated and real root-mean square errors (RMSE) for 
optimal and suboptimal (two-stage) algorithms (Nosov and 
Stepanov 2018). 


3 Accuracy Analysis Example 


Let us consider example of the gravity-aided navigation 
problem and compare the unconditional RMS errors for opti- 
mal and suboptimal (two-stage) algorithms. We have to spec- 
ify stochastic models for random process which describes the 
GA values and measurement errors. For the GA profile g(t) 
along a rectilinear trajectory we use the Jordan model in the 
form of a stationary process (Jordan 1972). Its correlation 
function can be written as follows: 


E co’) e—a 
5 ; 


where p > 0. The corresponding shaping filter is written as 
(Koshaev and Stepanov 2010; Peshekhonov and Stepanov 
2017): 


Kj (p) = o (: + ap (4) 


xı = —Bx, + X2, 

X2 = —Bx2 + x3, 

x3 = —Bx3 + qgwz, 
g = —Blx, + X2. 


(5) 


In Eqs. (4) and (5) p is a distance along the trajectory; a— 
inverse value of the correlation window; B = V oag/ap/ V2oz; 
V is the vessel speed; o93/a, is the parameter defining 
GA spatial variability; ggwz is forcing white noise with 
power-spectrum density (PSD) q% = 10°05; ož is the 


GA variance; and ¢ = (v5 = 1) / 5 is the dimensionless 
coefficient. 

To describe the errors of GA measurements on a sea ves- 
sel, we consider vertical movement due to heaving and white- 
noise error. The model describing vertical accelerations a, 
can be represented in the form: 


x4 = Xs, 
x5 = X6, 
ee oe ee: (6) 
X6 = —03X4 — 02X5 — D1 X6 + Wv, 
ay = X6, 
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Fig. 1 Example of the GA measurements and results of preliminary processing 


where b3 = (A? + w?)y3 b2 =A? + p? + 2uy; bi =2u +y; 
C = dyV2b3 (bib — b3) /bi, wy is white noise of unit 
PSD; o, is the RMS value of vertical movement x4; A is the 
prevailing heaving frequency; y is the coefficient of heaving 
irregularity; and y is the dimensionless coefficient. 

In this case, the gravimeter measurements can be written 
as 

y = Eta + gr = —PExi + x2 + X6 4+ Ver, (7) 
where vg, is the white-noise measurement error with a known 
PSD. 

Typical parameters for marine gravimetric survey were 
chosen for modeling the problem according to the discrete 
representation of Eqs. (5)-(7): the RMS value of the GA 
derivative was set to be 3 mGal/km, the period of vertical 
displacements was 7 s, and their RMS value was 0.2 m. 
Gravimeter measurements in the form Eq. (7) were simulated 
on a fixed section of several 30 km length GA profiles 
with coordinates (5,000 — 20,000) m. With spatial interval 
of 1 m it gives N = 15,000 raw measurements presented 
in Fig. la. Note that although the RMS value of the sea 
vessel vertical displacements is only 0.2 m, the RMS error 
of field measurements exceeds 14,700 mGal. It is obvious 
that against the background of such errors, it is impossible to 
estimate of the GA without processing. Figure 1b shows the 
results of preliminary processing in filtering and smoothing 
modes, as well as 30 bounds corresponding to the estimates. 
Values were obtained using the Kalman filter and the Rauch- 
Tung-—Striebel smoother based on models Eqs. (5)-(7). 


Using the GA estimates presented in Fig. 1b we can 
simplify the map-aiding by replacing model Eq. (6) in 
nonlinear part of the algorithm with white noise error model. 
It is feasible by subsampling the estimates, corrupted by 
correlated pre-processing errors. 

To select a subsampling interval, we use the correlation 
functions of pre-processing errors. The choice of subsam- 
pling interval based on value of 0.3 for the normalized 
correlation function. Under this condition, the spatial inter- 
vals for the pre-processed measurements were approximately 
1,200 and 800 m for the filtering and smoothing modes, 
respectively. The variance of the residual discrete white noise 
error was selected corresponding to the variance of the GA 
estimation error calculated during the pre-processing. 

The raw and pre-processed measurements, represented 
by Eqs. (2) and (3) were used in optimal and suboptimal 
algorithms, respectively, to simulate the solution of the 
gravity-aided navigation problem. In the optimal algorithm, 
a four-dimensional state vector including the A and vector 
Eq. (6) was estimated. The estimate was calculated as the 
mathematical expectation of the conditional probability den- 
sity function (p.d.f). In the same time suboptimal algorithms 
estimated only the A using pre-processed measurements. 
For calculations in both cases, we used the point-mass 
method with the number of nodes L = 3000, which was 
supplemented with the Rao-Blackwellization procedure for 
the optimal algorithm (Stepanov and Toropov 2015). The a 
priori RMS position error was set at 700 m. 

Figure 2 shows the examples of p.d.f. graphs for each 
algorithm. Orange lines indicate the true error of the NS, 
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Fig. 2 Examples of p.d.f. for optimal and suboptimal algorithms 


Table 1 Real RMSEs and calculation time for optimal and two-stage 


algorithms 
Number of measurements 
Algorithm Calc. time, s 
Optimal 302 
with filtering 
Suboptimal | 18 (800 m subsampling) 70 3.2 
with 


smoothing 


blue ones indicate its estimates obtained in the gravity-aiding 
algorithm, green ones indicate graphs of a posteriori p.d.f. 
depending on the distance covered. 

Graphs below show calculated and real unconditional 
RMSEs for the algorithms under study. They were calculated 
using 250 Monte Carlo simulations. Table | contains results 
including the values of real RMSE obtained at the end of the 
observation for n measurements and mean calculation time 
on the reference computer. 


Suboptimal with filtering 
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From the Fig. 3 and table above, it is obvious that the two- 
stage suboptimal algorithm with preliminary smoothing is 
comparable to the optimal algorithm in accuracy: when it was 
applied, the unconditional real RMSE of these algorithms 
were in 50-70 m range. In addition, calculated RMSE pro- 
vided by the suboptimal algorithm with smoothing is close 
to the real values. 

Algorithm with preliminary filtering of measurements 
underperforms both in accuracy and identity of real and cal- 
culated accuracy characteristics. Besides the smaller number 
of measurements used and greater variance of the estimation 
error, this can be explained by the presence of a phase delay 
in the field estimates produced by the filter. Although this 
algorithm is the fastest, low accuracy precludes its use. 

As for the amount of calculations, based on time aver- 
aging of 250 runs of algorithms on the test computer, we 
determined that it took 30 s to process all measurements in 
the optimal algorithm and an order of magnitude less, that is, 
3 s, for the two-stage algorithm which includes preliminary 
smoothing. 
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Fig. 3 Real and calculated unconditional RMSEs of the optimal and suboptimal algorithms 
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4 Conclusion 


The effect of preliminary gravity measurement processing on 
the accuracy of gravity-aided navigation has been analysed. 
It is based on comparison unconditional real RMS esti- 
mation errors for the two-stage suboptimal algorithms that 
use the measurement pre-processing with potential accuracy 
achieved by optimal Bayesian algorithm. The preliminary 
processing of the measurements is implemented in the fil- 
tering and smoothing modes. 

The example of gravity-aided navigation problem has 
been considered. It has been shown that the accuracy of the 
suboptimal algorithm with preliminary smoothing is close 
to the potential one, and the amount of calculations has 
been reduced by an order of magnitude in comparison with 
optimal algorithm. 

In the further research, it is supposed to consider map 
errors model and generalize the results to the case of update 
two-dimensional position of the vehicle. 
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Abstract 


The work considers the results of filtering and smoothing of the gravity disturbance 
vector horizontal components and focuses on the sensitivity of these results to the model 
parameters in the case when the inertial-geodesic method is applied in the framework of a 


marine survey on a sea vessel. 
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1 Introduction 


The complexity of the Earth’s surface and its interior is the 
reason why the direction of the gravity vector K4 (vertical 
line) does not coincide with the direction of the normal 
gravity vector T at points on the Earth’s physical surface. 
This difference is characterized by the gravity disturbance 
vector (GDV) (Peshekhonov et al. 2017; Torge 2001). The 
problem of determining all three components of this vector 
is a problem of vector gravimetry. Note that the vertical 
component of the GDV coincides with the free-air gravity 
anomaly up to the accuracy of correction for the height 
anomaly (Jekeli 1999, 2016). The horizontal components of 
the GDV are deviations of the vertical, which are determined 
by two angles of deviation: in the planes of the Meridian £ 
and the first vertical n (Peshekhonov et al. 2017; Torge 2001). 

Among the methods used to determine the GDV 
horizontal components on a moving base, including marine 
vessels, the gravimetric method, based on the measurement 
of gravity anomalies, and the astronomical-geodetic 
method, based on the comparison of astronomical and 
geodetic coordinates, and its modifications are particularly 
widespread (Koneshov et al. 2016; Peshekhonov et al. 
1995; Hirt et al. 2010). Also well-known are the method 
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of gravitational gradiometry, based on the measurement 
of the second derivatives of the geopotential, and the 
method of satellite or aircraft altimetry, which makes use of 
measurements of trajectory altitude (Koneshov et al. 2016). 
We cannot but mention methods based on the construction of 
global models of the Earth’s gravity field using the data 
on perturbations of satellite orbits (satellite-to-satellite 
tracking and other satellite missions). Combinations of 
these methods (for example, the astronomical-gravimetric 
method) are used (Peshekhonov et al. 2017; Koneshov et 
al. 2016). The inertial-geodetic method, which is one of the 
modifications of the astronomical-geodetic method, is used 
in the framework of the vector gravimetry problem solution. 
It is based on a comparison of astronomical latitude and 
longitude ọ, à, generated by the inertial navigation system 
(INS), and geodetic latitude and longitude B, L, obtained 
from the GNSS data (Koneshov et al. 2016; Emel’ yantsev et 
al. 2015). Thus, the horizontal components of the GDV are 
defined as: 


a es a) 
(à — L) = ņ/cos B. 
In contrast to the conventional astronomical-geodetic 
method, in the inertial-geodetic INS method, the INS gen- 
erates both astronomical coordinates and their derivatives, 
which makes it possible to use complementary velocity 
measurements (Koneshov et al. 2016; Dmitriev 1991). 
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The idea of using precision INS to determine horizontal 
components is based on the use of the relation between 
the errors of the INS output navigation parameters and the 
GDV components (Dmitriev 1997; Schultz and Winokur 
1969; Li and Jekeli 2008). The analysis of this error, which 
is actually a methodical one, opens up the possibility, in 
principle, of determining the GDV horizontal components 
directly relative to the reference value at a reference point 
by solving the estimation problem using differences in the 
measurements of INS and GNSS. 

The effectiveness of solving such a problem in real time 
(in filtering mode) was analyzed earlier in (Anuchin 1992; 
Nesenyuk et al. 1980; Staroseltsev and Yashnikova 2016). 
At the same time, experience shows that application of the 
smoothing mode, which does not require obtaining estimates 
of the desired parameters in real time, significantly increases 
the accuracy of the problem solution. Note that the feature 
of the smoothing problem is that when obtaining estimates 
at a current time moment, we can use both past and future 
measurements (Stepanov and Koshaev 2010; Loparev and 
Yashnikova 2012). The specifics of the inertial-geodetic 
method make it possible to implement optimal Kalman fil- 
tering and smoothing algorithms. However, such algorithms 
are based on the error models of the GNSS, INS and its 
sensors, as well as GDV horizontal components themselves 
in the form of random processes. To derive an adequate 
model of GDV horizontal components, it is necessary to 
have a priori information about the nature of the field in 
the survey area, for example, such parameters as the GDV 
variance and the correlation radius. Incorrect assignment 
of parameters can significantly reduce the accuracy of the 
resulting estimate. The analysis of the algorithm sensitivity 
is aimed to study the decrease in the estimation accuracy in 
the case that stochastic models in an algorithm are defined 
incorrectly (Stepanov and Koshaev 2003, 2011). Sensitivity 
of algorithms for determining the vertical component of 
the GDV was discussed in (Stepanov and Motorin 2019; 
Stepanov et al. 2015a, b). 

The paper considers the results of filtering and smoothing 
of the GDV horizontal components and focuses on the 
sensitivity of these results to the model parameters in the case 
when the inertial-geodetic method is applied in the frame- 
work of a marine survey on a sea vessel. The studies involve 
direct calculation of the covariance matrices of estimation 
errors (without using the Monte Carlo method). 


2 Description of the Error Model 


To solve the problems of filtering and smoothing in this 
study, we used the model of INS errors described in (Gusin- 
sky et al. 1996, 1997). The state vector of this model 
includes errors in constructing the local vertical, errors in 
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generation of the East and North velocity components, errors 
of the sensing elements, and the model of GDV horizontal 
components. Loosely coupled INS-GNSS system is assumed 
with coordinate, velocity, orientation angle, and sensor errors 
feedback for the filtering mode during survey. The fixed 
interval smoothing for every tack was fully done after survey. 

The approach to creation of an integrated system to 
measure GDV horizontal components implies very high 
requirements imposed on the instrumental errors of the INS. 
Gyroscope drift and measurement errors of linear accelera- 
tion must not exceed 0.0001 deg./h and 1 arcsec, respectively 
(Emel’yantsev et al. 2015; Nesenyuk et al. 1980). Therefore, 
the coefficients of the gyroscope error model were described 
by first-order Markov processes with the standard deviation 
(SD) at a level of 3 - 1075°/h and a correlation interval of 
40 h; nonwhite-noise errors of the accelerometers were set 
by the stationary first-order Markov processes with an SD 
at a level of 1 arcsec and a correlation interval of 10 h; the 
additive white noise of the accelerometers was also set at a 
level of 1 arcsec/,/Hz. These processes define stability of 
inertial sensors. The initial root mean square errors (RMSE) 
of the sensors were set equal to the steady-state values of the 
corresponding Markov process SD. The RMSE of the INS 
initial alignment were set at a level of 1 arcsec for tilt angles 
and at a level of 0.1 m/s for the horizontal components of the 
velocity. 

As noted above, for the solution of the estimation problem 
by the inertial-geodetic method in the formulation consid- 
ered here, models of the GDV horizontal components must 
be given in the form of random processes describing the 
variability of the components along the vehicle trajectory 
(Staroseltsev and Yashnikova 2016; Nash and Jordan 1978). 
The expressions of spectral densities for the longitudinal 
S-i(@) and transverse S,7(@) horizontal components of the 
GDV in the rectilinear motion of the vehicle with a constant 
velocity V were assumed to be similar to those given in 
(Nesenyuk et al. 1980): 


4a,02w2 2a,02 V 
Se = Z $ Se = ————— ; e= ai 
L (w) (w? + a2} ate) w +a? Oe 93d 
(2) 


where Os, Œs, d are SD, damping coefficient, and inverse 
correlation radius, respectively. 

For the solution of the estimation problem, the model 
included measurements from the INS and the GNSS receiver, 
which are parts of the integrated system. It was assumed 
that the position measurements of the GNSS had white noise 
at a level of 3 cm/,/Hz and a component represented by 
a stationary first-order Markov process with a SD of 3 cm 
and a correlation interval of 1 h. Velocity measurements 
of the GNSS had only a white-noise error at a level of 
3 cm/s/ V/Hz. 
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3 Estimating the Accuracy 
of Determining the GDV Horizontal 
Components: The Results 


To simulate the solution of the smoothing problem, we 
modified software in the Matlab package (Stepanov and 
Koshaev 2003, 2011), which allowed us to calculate the 
RMSE of filtering and smoothing estimates and display their 
plots. Error-tolerant algorithms were used to calculate error 
covariance matrices based on UD decomposition. Owing to 
this, all the necessary calculations were performed without 
additional computational errors. 

In the simulation, it was assumed that the vessel was mov- 
ing North along the meridian at a speed of 10 knots; the initial 
inertial longitude was 30° E; the initial latitude was 60° N. 
The solution interval was 24 h, and the discreteness — 10 s. 

We set different characteristics of the gravitational field 
variability: from weakly disturbed with a of 5 arcsec and 
a correlation radius of 20 nmi to the highly disturbed with 
an RMSE of 20 arcsec and a correlation radius of 5 nmi. 
The steady-state values of the RMSE in the estimation of the 
GDV horizontal components obtained in the simulation for 
filtering and smoothing algorithms are shown in Table 1. 

Using the smoothing mode to calculate the GDV hori- 
zontal components provides a significant gain in accuracy 
compared to the filtering mode. As may be seen from the 
table, for a weakly disturbed field (SD of 5 arcsec) only 
filtering algorithms can be used because there is practically 
no gain from smoothing. This fact is easy to explain: it is 
well known that the estimation accuracy of constant values 
(random bias) in the filtering and smoothing modes are the 
same. As for a strongly disturbed field (SD of 20 arcsec), 
the use of smoothing algorithms can give almost a two-fold 
increase in the accuracy of the estimation problem solution. 
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4 Estimating the Sensitivity 
of Algorithms for Determining 
the GDV Horizontal Components 
to the Model Parameters: The Results 


For sensitivity analysis, we used the same parameters of the 
GDV model as in Table 1. When calculating the sensitivity 
for each combination of real (true) and estimated values of 
a model parameter, we determined three different RMSE in 
the estimation of the GDV horizontal components. 

The first one is the RMSE of the error in estimating the 
optimal algorithm, which corresponds to the estimate in the 
case that the model specified in the algorithm corresponds 
to the real one. This RMSE characterizes the maximum 
achievable estimation accuracy of the problem under given 
conditions. 

The second one is the real RMSE that correspond to 
the real estimation accuracy of the algorithm tuned to an 
incorrect model of a disturbed field. Real RMSE calculated 
using method from (Stepanov and Koshaev 2003, 2011). The 
method is based on the solution of the covariance equation 
for the augmented vector, which includes both optimal and 
suboptimal estimates of state vectors corresponding to cor- 
rect and incorrect models. This RMSE will always be higher 
than the RMSE for the optimal algorithm. 

The third RMSE are those calculated using the data from 
the covariance equation of the algorithms that stand for the 
characteristics of the accuracy of the obtained solutions. If 
for different values of the parameters of the real models in a 
certain range the calculated RMSE is not lower than the real 
one, it may be stated that the algorithm has a guaranteeing 
property. We assume suitable to provide such RMSE along 
with others because it is one of the ways to judge about 
accuracy in real survey. 


Table 1 The steady-state values of the RMSE in the estimation of the GDV horizontal components obtained in the simulation for filtering and 


smoothing algorithms 


Field SD (arcsec) Correlation radius (nmi) GDV component Filtering RMSE (arcsec) | Smoothing RMSE (arcsec) 
Aé — North (longitudinal) 0.8 0.7 
5 20 An — East (transversal) 1.2 1.1 
A£ — North (longitudinal) 1.2 0.8 
10 15 An — East (transversal) 1.4 1.2 
A£ — North (longitudinal) 1.7 1.0 
15 10 An — East (transversal) 1.6 1.2 
Aé — North (longitudinal) 2.6 1.3 
20 5 An — East (transversal) 2.2 1.4 
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Fig. 1 RMSE of the GDV horizontal components for the case of motion in a strongly disturbed field (field SD — 20 arcsec, correlation radius — 
5 nmi) with the filter and smoother tuned to a weakly disturbed field (field SD — 5 arcsec, correlation radius — 20 nmi) 


An example of simulation for optimistic tuning of the 
algorithm, that is, such that the algorithm is tuned to a less 
variable field is shown in Fig. 1. It is a case of motion 
in a strongly disturbed field (SD of the field is 20 arcsec, 
correlation radius — 5 nmi) with a setting to a weakly 
disturbed field (SD of the field — 5 arcsec, correlation radius — 
20 nmi). 

Analyzing the curves in Fig. 1, we can draw the following 
conclusions: for the motion in a highly disturbed field with 
the filter, which use incorrect weakly disturbed model. The 
estimation accuracy will be approximately 1.5 times worse 
than the optimal one and 2-3 times worse than the calculated 


one, which will be shown by the covariance equation of 
the suboptimal filter. These conclusions are also valid for 
the case of smoothing algorithms. However, due to the fact 
that the smoothing mode is more accurate than the filtering 
mode, the difference in the absolute values of the RMSE 
in estimating the GDV orizontal components between the 
suboptimal and optimal algorithms is not so significant here. 
As can be seen from Fig. 1, for the preset parameters, the 
guaranteeing property of the filter does not hold either. 

The situation that arises is not satisfactory for estimation 
of the GDV horizontal components, which is why an attempt 
was made to find such a filter setting that would ensure min- 
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Fig. 2 RMSE of GDV horizontal components for the case of motion in a strongly disturbed field (field SD — 20 arcsec correlation radius — 5 nmi) 
with the filter and smoother tuned to a field with SD of 15 arcsec and correlation radius of 10 nmi 


imal loss to the optimal algorithm both in weakly disturbed 
and strongly disturbed fields. Figure 2 presents similar plots 
for the motion in a strongly disturbed field with a SD of 
20 arcsec and a correlation radius of 5 nmi with the filter 
setting for the field SD of 15 arcsec and a correlation radius 
of 10 nmi. 

With this setting, the suboptimal filter insignificantly loses 
to the optimal one; however, the guaranteeing property of 


estimation does not hold either. The results of the simulation 
of motion in a weakly disturbed field (the field SD — 5 arcsec, 
the correlation radius — 20 nmi) with the same setting are 
presented in Fig. 3. 

It is obvious that in this case, too, the real RMSE of 
the error in estimating the GDV horizontal components is 
close to the optimal one. However, the calculated RMSE 
significantly differs from both the optimal and real RMSE. 
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Fig. 3 RMSE of GDV horizontal components for the case of motion in a weakly disturbed field (field SD — 5 arcsec, correlation radius — 20 nmi) 
with the filter and smoother tuned to a field with SD of 15 arcsec and the correlation radius of 10 nmi 


5 Conclusion 


The accuracy of the problem solution in determining GDV 
horizontal components by the inertial-geodetic method on a 
marine moving vessel in the filtering and smoothing modes 
has been studied. It is shown that the use of the smoothing 
mode for estimation of GDV horizontal components leads 
to a significant gain in accuracy in a strongly disturbed field 
as compared to the filtering mode. However, for a weakly 
disturbed field, there is practically no gain from smoothing. 
Sensitivity of the geodetic method algorithms for estimat- 
ing GDV horizontal components to the accuracy of setting 
the model of these components was analyzed. The analysis 
has shown that with optimistic filter tuning (tuning to a field 
that is less variable than in reality), suboptimal algorithms 
can significantly lose in accuracy to the optimal algorithm, 
and they do not provide guaranteeing estimation. At the same 
time, within the known limits of changes in the correlation 
radiuses and the RMSE of the field in the survey area, 
it is possible to choose the model parameters so that the 
suboptimal algorithm will not differ noticeably in accuracy 


from the optimal one. However, it is difficult to obtain an 
adequate calculated characteristic of accuracy. The solution 
to this problem might be in the development of adaptive 
algorithms. 
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